genetic algorithm assignment problem

An official website of the United States government

The .gov means it’s official. Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

The site is secure. The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

• Publications
• Account settings
• Advanced Search
• Journal List
• Entropy (Basel)

A Hybrid Genetic-Hierarchical Algorithm for the Quadratic Assignment Problem

Associated data.

Not applicable.

In this paper, we present a hybrid genetic-hierarchical algorithm for the solution of the quadratic assignment problem. The main distinguishing aspect of the proposed algorithm is that this is an innovative hybrid genetic algorithm with the original, hierarchical architecture. In particular, the genetic algorithm is combined with the so-called hierarchical (self-similar) iterated tabu search algorithm, which serves as a powerful local optimizer (local improvement algorithm) of the offspring solutions produced by the crossover operator of the genetic algorithm. The results of the conducted computational experiments demonstrate the promising performance and competitiveness of the proposed algorithm.

1. Introduction

The quadratic assignment problem (QAP) [ 1 , 2 , 3 , 4 , 5 , 6 ] is mathematically formulated as follows: Given two non-negative integer matrices A = ( a i j ) n × n , B = ( b k l ) n × n , and the set Π n of all possible permutations of the integers from 1 to n , find a permutation p = ( p ( 1 ) , p ( 2 ) , … , p ( n ) ) ∈ Π n that minimizes the following objective function:

One of the examples of the applications of the quadratic assignment problem is the placement of electronic components on printed circuit boards [ 7 , 8 ]. In this context, the entries of the matrix A are associated with the numbers of the electrical connections between the pairs of components. The entries of the matrix B correspond to the distances between the fixed positions on the board. The permutation p = ( p ( 1 ) ,   p ( 2 ) ,   … ,   p ( n ) ) can be interpreted as a separate configuration for the arrangement of components in the positions. The element p ( i ) in this case indicates the number of the position to which the i -th component is assigned. In this way, z (or more precisely, z / 2 ) can be understood as the total (weighted) length of the connections between the components, when all n components are placed into the corresponding n positions.

The other important areas of applications of the QAP are as follows: assigning runners in relay teams [ 9 ], clustering [ 10 ], computer/telephone keyboard design [ 11 , 12 , 13 ], planning of airport terminals [ 14 ], facility location [ 15 ], formation of chemical molecular compounds [ 16 ], formation of grey patterns [ 17 ], index assignment [ 18 ], microarray layout [ 19 ], numerical analysis [ 20 ], office assignment and planning of buildings [ 21 , 22 ], seed orchard design [ 23 ], turbine balancing [ 24 ], website structure design [ 25 ]. More examples of the practical applications of the QAP can be found in [ 4 , 26 ].

The quadratic assignment problem is also a complicated theoretical-mathematical problem. It is proved that the QAP belongs to the class of the NP-hard optimization problems [ 27 ]. The QAP can be solved exactly if the problem size ( n ) is quite small ( n < 30 ) [ 28 , 29 , 30 , 31 , 32 , 33 , 34 , 35 ]; although some special case QAP examples of larger sizes ( n = 36 [ 36 ], n = 64 [ 10 ]) have also been exactly solved. For this reason, heuristic optimization algorithms are widely used. Although these algorithms do not guarantee the optimality of the obtained solutions, they allow for a sufficiently high quality (near-optimal) solutions within a reasonable computation time [ 37 ].

Classical single-solution local search (LS) and related algorithms were intensively used for the QAP in the early period of application of heuristic algorithms (1960–1980) [ 38 , 39 , 40 , 41 ]. Later, improved local search algorithms have been employed [ 42 , 43 , 44 , 45 ]. The so-called breakout local search has been empirically proven to be quite efficient [ 46 , 47 ].

Simulated annealing (SA)-based algorithms usually provide better quality results, compared to pure, deterministic local search algorithms. This applies to both the early variants of SA algorithms [ 48 , 49 , 50 ] and improved SA algorithm modifications [ 51 , 52 , 53 ].

Even more performance was achieved by adopting the tabu search (TS) methodology-based algorithms. The fast-running tabu search algorithm developed by Taillard [ 54 ] in the early 1990s is still considered as one of the most successful single-solution-based heuristic algorithms for the QAP. Since then, a number of improved variations of TS algorithms have been proposed: reactive tabu search [ 55 ], concentric tabu search [ 56 ], enhanced tabu search [ 57 ], tabu search with hardware acceleration [ 58 ], self-controlling tabu search [ 59 ], repeated iterated tabu search [ 60 , 61 ], parallel tabu search [ 62 ], and other variants [ 58 , 63 ]. The performance of LS and TS algorithms can be increased by extending these algorithms to their ameliorated “siblings”, namely, the iterated LS (ILS) [ 64 , 65 ] and iterated TS (ITS) algorithms [ 66 ]. Iterated search algorithms have some similarities with multistart methods [ 63 , 67 , 68 , 69 ] as well as greedy adaptive search procedures (GRASPs) [ 70 ].

Population-based/evolutionary algorithms constitute another important class of efficient heuristic algorithms for the QAP. The advantage of this class of algorithms is that these algorithms operate with the sets of solutions instead of single solutions and this property is of prime importance when it comes to the solution of the QAP and related problems. In particular, it is found that, namely, the genetic algorithms (GA) seem to be very likely among the most powerful heuristic algorithms for solving the QAP, among them: greedy genetic algorithm [ 71 ], genetic-local search algorithm [ 72 , 73 , 74 ], genetic algorithm using cohesive crossover [ 75 ], improved genetic algorithm [ 76 ], parallel genetic algorithm [ 77 ], memetic algorithm [ 78 ], genetic algorithm on graphics processing units [ 79 ], quantum genetic algorithm [ 80 ], and other GA modifications [ 81 , 82 , 83 , 84 , 85 , 86 , 87 , 88 , 89 ]. Note that the population-based algorithms are usually hybridized with the single-solution-based algorithms (local search, tabu search, iterated local/tabu search, GRASP). Overall, a significant part of the algorithms for the QAP are, in essence, hybrid algorithms [ 71 , 73 , 75 , 76 , 78 , 79 , 82 , 83 , 88 , 90 , 91 , 92 , 93 , 94 , 95 , 96 , 97 , 98 , 99 , 100 , 101 , 102 , 103 , 104 ]. It is the hybrid genetic-iterated tabu and genetic-breakout local search algorithms [ 76 , 78 ] that allowed to achieve the most promising results.

Swarm intelligence algorithms simulate collective intelligent behaviour of physical/biological entities (agents) (like particles (particle swarm optimization algorithms [ 105 , 106 ]), ants (ant colony optimization algorithms [ 107 ]), bees (artificial bee colony algorithms [ 108 , 109 ]). Finally, the algorithms inspired from real-world phenomena (including those using metaphors) are also applicable to the QAP [ 90 , 96 , 98 , 110 , 111 , 112 , 113 , 114 , 115 , 116 , 117 , 118 ]. For more extensive surveys and literature reviews on the QAP, the reader is referred to [ 4 , 119 ].

The main contribution of this article is that it presents an innovative hierarchicity-based genetic algorithm which is hybridized with a multi-level hierarchical iterated tabu search (HITS) algorithm serving as a powerful optimizer of the offspring solutions. The basic idea of HITS is, in turn, based on the multiple (re)use of the iterated tabu search (ITS) and, simultaneously, moving through many different locally optimal solutions. The other important novelty is that the original crossover and mutation operators are introduced. The crossover operator distinguishes for its universality and, at the same time, versatility and flexibility; while the mutation operation is integrated within the HITS algorithm and is based on combined deterministic and probabilistic (controlled random) moves between solutions during the tabu search process. Also, we have employed the modified population replacement rule. Finally, we have incorporated the population restart mechanism to avoid search stagnation. All these new features have led to the development of high-performance genetic algorithm with excellent results.

The remainder of the paper is structured as follows: In Section 2 , some basic definitions are given. Then, in Section 3 , the detailed description of the genetic-hierarchical algorithm and its constituent parts is provided. In Section 4 , the results of the computational experiments with the proposed algorithm are presented. The paper is completed with concluding remarks.

2. Basic Definitions

At the beginning, we provide some preliminary (basic) formal definitions.

Suppose that p ( u ) ( u = 1 , … ,   n ) and p ( v ) ( v = 1 , … ,   n , u ≠ v ) are two elements of the permutation p . Then p u , v is defined in the following way:

This means that the permutation p u , v is obtained from the permutation p by interchanging exactly the elements p ( u ) and p ( v ) in the permutation p . The formal operator (move, or transition operator) ϕ ( p , u , v ) : Π n × N × N → Π n swaps the u th and v th elements in the given permutation such that p u , v = ϕ ( p , u , v ) . Note that the following equations hold: p u , u = p , p u , v = p v , u , ( p u , v ) u , v = p .

The difference in the objective function ( z ) values when the permutation elements p ( u ) and p ( v ) are interchanged is calculated according to the following formula:

The difference in the objective function values can be calculated more faster—under condition that the previously calculated differences ( Δ ( p i , j ,   p ) ( i ,   j = 1 , … , n )) are memorized (stored in a matrix Ξ ). The difference value is calculated using O ( 1 ) operations [ 54 , 120 ]:

After the interchange of the elements p ( u ) and p ( v ) has been performed, new differences Ξ ′ ( i , j ) ( i ,   j = 1 , … , n , i ≠ u , i ≠ v , j ≠ u , j ≠ v ) are calculated according this equation:

If i = u or i = v or j = u or j = v , then the Formula (3) should be applied. So, all the differences are calculated using only O ( n 2 ) operations. Still, the initial calculation of the values of the matrix Ξ requires O ( n 3 ) operations (but only once before starting the optimization algorithm).

If the matrix A and/or matrix B are symmetric, then Formula (3) becomes simpler. Assume that the matrix B is symmetric. Then, the (asymmetric) matrix A can be transformed to symmetric matrix A ′ . Thus, we get the following formula:

here, a ′ u k = a u k + a k u , u = 1 , … , n , v = 1 , … , n , u ≠ v . Analogously, Formula (5) turns to the formula:

If i = u ( v ) or j = u ( v ) , Equation (6) must be applied.

In addition to this, suppose that we dispose of 3-dimensional matrices A ″ = ( a ″ u v k ) n × n × n and B ″ = ( b ″ l r t ) n × n × n . Also, let a ″ u v k = a ′ u k − a ′ v k , b ″ l r t = b l t − b r t , l = p ( j ) , r = p ( i ) , t = p ( k ) . Then, we can apply the following formulas for calculation of the differences in the objective function values:

Using the matrices A ″ and B ″ allows to save up to 20 % of computation (CPU) time [ 66 ]. The distance (Hamming distance) between two permutations p and p ′ is defined as:

The following equations hold: δ ( p ,   p ) = 0 , δ ( p ,   p ′ ) ≠ 1 , 0 ≤ δ ( p ,   p ′ ) ≤ n , δ ( p ,   p ′ ) = δ ( p ′ ,   p ) , δ ( p ,   p u , v ) = 2 for any u , v ( u ≠ v ). In the case of disposing of k different numbers u 1 ,   u 2 ,   … ,   u k , this equation holds: δ ( p , ( ( ( p u 1 ,   u 2 ) u 2 ,   u 3 ) … ) u k − 1 ,   u k ) = k .

The neighbourhood function Θ : Π n → 2 Π n assigns for each p ∈ Π n its neighbourhood (the set of neighbouring solutions) Θ ( p ) ⊆ Π n . The 2-exchange neighbourhood function Θ 2 is defined in the following way:

where δ ( p ,   p ′ ) is the Hamming distance between solutions p and p ′ . The neighbouring solution p ′ ∈ Θ 2 ( p ) can be obtained from the current solution p by using the operator ϕ ( p ,   ⋅ ,   ⋅ ) . The computational complexity of exploration of the neighbourhood Θ 2 is proportional to O ( n 2 ) .

Solution p l o c _ o p t ∈ Π n is said to be locally optimal with respect to the neighbourhood Θ if for every solution p ′ from the neighbourhood Θ ( p l o c _ o p t ) the following inequality holds: z ( p l o c _ o p t ) ≤ z ( p ′ ) .

3. Hybrid Genetic-Hierarchical Algorithm for the Quadratic Assignment Problem

Our proposed genetic algorithm (for more thorough information on the genetic algorithms, we refer the reader to [ 121 ]) is based on the hybrid genetic algorithm framework, where explorative (global) search is in tandem with the exploitative (local) search. The most important feature of our genetic algorithm is that it is hybridized with the so-called hierarchical (self-similar) iterated tabu search (HITS) algorithm (see Section 3.4 ).

The permutation elements p ( 1 ) ,   p ( 2 ) ,   … ,   p ( n ) are directly linked to the individuals’ chromosomes—such that the chromosomes’ genes correspond to the single elements p ( 1 ) ,   p ( 2 ) ,   … ,   p ( n ) of the solution p . No encoding is needed. The fitness of the individual is associated with the corresponding objective function value of the given solution, z ( p ) .

The following are the main essential components (parts) of our genetic-hierarchical algorithm: (1) initial population construction; (2) selection of parents for crossover procedure; (3) crossover procedure; (4) local improvement of the offspring; (5) population replacement; (6) population restart. The top-level pseudo-code of the genetic-hierarchical algorithm is presented in Algorithm 1 (Notes: (1) The subroutine GetBestMember returns the best solution of the given population. (2) The mutation process is integrated within the k -level hierarchical iterated tabu search algorithm. The mutation process depends on the mutation variant parameter MutVar .).

3.1. Creation of Initial Population

There are two main population construction phases. In the first one, the pre-initial population is constructed and improved; in the second one, the culling of the improved population is performed. So, firstly, P S ′ = P S × C 1 members of the pre-initial population P are created using the version of the GRASP algorithm [ 70 ] implemented by the authors. P S denotes the population size, and C 1 ( C 1 ≥ 1 ) is the user-defined parameter and is to regulate the size of the pre-initial population.

There are several options of the population construction in the first phase controlled by the parameter I n i t P o p V a r . If I n i t P o p V a r = 1 , then every generated solution is improved by the hierarchical iterated tabu search algorithm. There are few conditions. If the improved solution ( p ✩ ) is better than the best so far found solution in the population P , then the improved solution replaces the best found solution. Otherwise, it is tested if the minimum mutual distance between the improved solution ( p ✩ ) and the existing population members ( min p ∈ P { δ ( p ✩ ,   p ) } ) is greater than or equal to the predefined distance threshold, D T . If this is the case, the improved solution is added to the population P . Otherwise, the improved solution is disregarded and simply a random solution is added instead. (Remind that the distance between solutions is calculated using Equation (10)). The distance threshold is obtained from the following equation: D T = max { 2 ,   ⌊ ε n ⌋ } , where ε denotes the distance threshold factor ( 0 < ε ≤ 1 ). This presented scheme is to ensure the high level of diversity of the population members and, at the same time, to enhance the searching ability of the genetic algorithm. To obtain better initial population, the HITS algorithm with the increased number of iterations is used during the initial population formation. This is similar to a compounded approach proposed in [ 122 ].

The second option ( I n i t P o p V a r = 2 ) is almost identical to the first one, except that the genetic algorithm itself (a de facto copy of the hybrid genetic-hierarchical algorithm) (instead of the HITS algorithm) is employed for the creation of the initial population. As an alternative option ( I n i t P o p V a r = 3 ) of the population improvement, two-level genetic-hierarchical algorithm (master-slave genetic algorithm) can be employed for the initial population improvement.

In the second phase—which is very simple— ( C 1 − 1 ) P S worst members of the pre-initial population are truncated and only P S best members survive for the subsequent generations.

3.2. Selection of Parents

The selection of parents is performed by using the parametrized rank-based selection rule [ 123 ]. In this case, the positions ( κ 1 , κ 2 ) of the parents within the sorted population are determined according to the following formulas: κ 1 = ⌊ ( ς 1 ) σ ⌋ , κ 2 = ⌊ ( ς 2 ) σ ⌋ , κ 1 ≠   κ 2 , where ς 1 , ς 2 are uniform (pseudo-)random numbers in the interval [ 1 ,   P S 1 σ ] , here P S denotes the population size, and σ is a real number from the interval [ 1 ,   2 ] (it is referred to as a selection factor). It is clear that the better the individual, the larger the selection probability.

3.3. Crossover Operators

Two-parent crossover is formally defined by using operator Ψ :   Π n × Π n → Π n such that:

where p ′ , p ″ , p ° denote parental solutions, and p ° is the offspring solution. (The child can coincide with one of the parents if, for example, the parents are very similar.) The crossover operator must ensure that the chromosome of the offspring necessarily inherits those genes that are common in both parent chromosomes, i.e., (also see Figure 1 ):

here, p ′ , p ″ , p ° refer to the parents and the offspring, respectively.

Graphical illustration of a crossover.

In our work, the crossover procedure takes place at every generation of the genetic-hierarchical algorithm, i.e., the crossover probability is equal to 1 . Several crossover operators were implemented and examined. Short descriptions of the crossover procedures are provided below (see also [ 124 , 125 ]).

3.3.1. One-Point Crossover—1PX

1PX is among the most popular genetic crossover operators. Very briefly, the basic idea is as follows. A single crossover point (position, or locus) is chosen in one of the two chromosomes. The position x can be determined by generating a uniformly distributed (pseudo-)random number within the interval [ 1 ,   n − 1 ] ( n is the chromosome length). The offspring is obtained by copying x genes from one parent, the rest of genes are copied from the opposite parent. If there are empty loci left, they are filled in randomly; in addition, the feasibility of the offspring must be preserved.

3.3.2. Two-Point Crossover—2PX

Two-point crossover works similarly to the one-point crossover, except that two crossover points x 1 and x 2 ( 1 ≤ x 1 < x 2 < n ) are used.

3.3.3. Uniform Crossover—UX

In this case, the common genes are copied to the offspring’s chromosome. Then, the unoccupied positions in the offspring’s chromosome are scanned form left to right and the empty loci are assigned the genes—one at a time—from one of the parents with probability 1 2 , i.e., p ° ( i ) = { p ′ ( i ) , ς < 1 2 p ″ ( i ) , otherwise ; here, ς is a (pseudo-)random number from the interval [ 0 ,   1 ] . The assigned gene must be unique.

3.3.4. Shuffle Crossover—SX

The shuffle crossover is obtained by randomly reordering the parents’ genes before applying the uniform crossover. The same rearrangement rule must be used for both parents. After the uniform crossover is finished, the same (initial) rearrangement rule is again applied.

3.3.5. Partially-Mapped Crossover—PMX

Partially-mapped crossover can be seen as a modified variant of the k -point (multi-point) crossover. The basic principle relies on the so-called mapping sections (the chromosome segments between mapping points). So, at first, the segments of the chromosome of one parent are moved to the offspring’s chromosome. The same is done for the other parent. At last, the empty loci (if any) are filled in in a random way.

3.3.6. Swap-Path Crossover (SPX)

3.3.6.1. swap-path crossover (basic version)—spx1.

The main distinguishing feature of SPX is that instead of transferring genes from parents to a child, the genes are, so to speak, rearranged in the chromosomes of the parents. Let ( p ′ ,   p ″ ) be a pair of parents. Then, the process starts from an arbitrary position and the positions are scanned from left to right. The process continues until a predefined number of swaps, s ( s < n ), have been performed. If, in the current position, the genes are the same for both parents, then one moves to the next position; otherwise, a pairwise interchange of genes of the parents’ chromosomes is accomplished. The interchange is performed in both parents. For example, if the current position is i and a = p ′ ( i ) , b = p ″ ( i ) , then there exists a position j such that b = p ′ ( j ) , a = p ″ ( j ) ; then, after a swap, p ′ ( i ) = b , p ″ ( i ) = a and p ′ ( j ) = a , p ″ ( j ) = b . Consequently, new chromosomes, say p ‴ , p ⁗ , are produced. In the next iteration, a pair ( p ‴ ,   p ⁗ ) is considered, and so on.

3.3.6.2. Swap-Path Crossover (Modified Version I)—SPX2

This modification is achieved when the best offspring (with respect to the fitness of the offspring) is retained in the course of gene interchanges.

3.3.6.3. Swap-Path Crossover (Modified Version II)—SPX3

The essential feature this crossover procedure is that the offspring fitness is dynamically evaluated: only the gene interchanges that improve the value of the objective function are accepted.

3.3.7. Cycle Crossover—CX

The cycle crossover is based on the pairwise gene interchanges. The key property of CX is the ability to produce the offspring without distortion of the genetic code; in the other words, CX enables to create the chromosome with no random (foreign) genes. The negative aspect of CX is that the offspring may genetically be very close to their predecessors.

3.3.8. Cohesive Crossover—COHX

The cohesive crossover was proposed by Z. Drezner [ 75 ] to efficiently recombine individuals’ genes by taking into account the problem data, in particular, the distances between objects’ locations. From several recombinations of genes, the recombination is selected that minimizes the objective function.

3.3.9. Multi-Parent Crossover—MPX

In the multi-parent crossover, several (or all) members of a population participate in creation of the offspring. More precisely, the i th position (locus) of the offspring’s chromosome p ° is assigned the value j with the probability P ( p ( i ) = j ) (under condition that the value j has not been utilized before).

The probability that p ( i ) = j ( P ( p ( i ) = j ) ) is calculated according to the formula: P ( p ( i ) = j ) = q i j ∑ j = 1 n q i j ; where q i j is an element of the matrix Q = ( q i j ) n × n ; here, q i j denotes the number of times that the i th locus takes the value j in the parental chromosomes. If there exist several values ( j 1 , j 2 , …) with the same probability, then one of them is chosen randomly.

3.3.10. Universal Crossover—UNIVX

The universal crossover (UNIVX) [ 124 ] distinguishes for its versatility and the possibility of flexible usage depending on the specific needs of the user. It is somewhat similar to what is known as a simulated binary crossover [ 126 ].

Our operator is based on the use of a random mask. There are four controlling parameters: χ 1 , χ 2 , χ 3 , χ 4 . The mask length is equal to χ 1 , where χ 1 is a (pseudo-)random number within the interval [ ε 1 ,   n ] , n is the length of the chromosome, ε 1 = ⌊ r × n ⌋ , r is the user’s parameter close to 1 , for example, 0.9 . The mask contains binary values 0 and 1 , where 1 signals that the corresponding gene of the first parent’s chromosome must be chosen and 0 is to indicate that the second parent’s gene needs to be replicated. The degree of randomness of the mask is controlled by the parameters χ 2 , χ 3 . The parameter χ 2 ( χ 2 ∈ [ ε 2 ,   ε 3 ] , 0 < ε 2 ≤ ε 3 < 1 ) dictates how many 0 ’s and 1 ’s are there in the mask: the higher the value of χ 2 , the bigger total number of 1 ’s, and vice versa. The juxtaposition of bits is regulated by the parameter χ 3 . The bit generation itself is performed by using a kind of “anytime” min-max sorting algorithm. That is, if the sorting algorithm is interrupted at some random moment, this results in a randomized (“quasi-sorted”) sequence of bits. The moment of interruption is defined by the number η , where η = χ 3 w , here χ 3 is a (pseudo-)random real number from the interval [ 0 ,   1 ] , and w denotes the maximum number of iterations required to fully sort all the bits. (As an example, if the bits “ 0000001111 ” are to be sorted in the descending order and the algorithm is stopped at χ 3 = 0.9 , then the random mask similar, for example, to “ 101100010 0” would be generated.) Having the mask generated, the decision is made as to about what genes have to be transmitted to the offspring’s chromosome. The index of the starting locus of the transferred genes, χ 4 , is generated randomly—in such a way that χ 4 ∈ [ 1 ,   n ] . Eventually, the empty loci (if any) are filled in randomly.

3.4. Offspring Improvement

3.4.1. hierarchical iterated tabu search algorithm.

Every created offspring is improved by using the hierarchical iterated tabu search algorithm, which inspires from the philosophy of iterated local search [ 127 ] and also the spirit of self-similarity—one of the fundamental properties of nature (see Figure 2 ). Basically, this means that the algorithm is (almost) exactly similar to the part of itself. In the other words, the main idea is the repeated, cyclical adoption (reuse) of the iterated tabu search algorithm, that is, the iterated tabu search can be reused multiple times. This idea is not very new, and some variants of hierarchical-like algorithms have been already investigated [ 128 , 129 , 130 , 131 , 132 , 133 , 134 ].

Visual conceptual interpretation of hierarchicity.

The paradigm of the hierarchicity based (self-similar) algorithm is as follows:

• (1) Set the number of levels, k ( 1 ≤ k ≤ k m a x ).
• (2) Generate an initial solution p .
• (3) Apply k ‒1-level algorithm to the solution p . Let p ✸ be the improved solution.
• (4) Memorize the best found solution.
• (5) Set p ☼ = p ✸ or select a solution p ☼ from the history of solutions.
• (6) Apply a perturbation procedure to the solution p ☼ . Let p ~ be the perturbed solution.
• (7) Set p = p ~ .
• (8) If the termination criterion is not satisfied, then go to Step 2; otherwise, stop the algorithm.

The k -level hierarchical iterated tabu search algorithm consists of three basic components: (1) invocation of the k –1-level hierarchical iterated tabu search algorithm to improve a given solution; (2) acceptance of the candidate (improved) solution for perturbation, i.e., mutation; (3) mutation of the accepted solution.

Perturbation (mutation) is applied to a chosen optimized solution that is selected by the defined candidate acceptance rule (see Section 3.4.3 ). The mutated solution serves as an input for the self-contained TS procedure. The TS procedure returns an improved solution, and so on. The overall process continues until a pre-defined number of iterations have been performed (see Algorithm 2 (Note. The iterated tabu search procedure (see Algorithm 3) is in the role of the 0-level hierarchical iterated tabu search algorithm.)). The best found solution is regarded as the resulting solution of HITS.

The 0-level HITS algorithm is in fact simply iterated tabu search algorithm (for more details on the ITS algorithm, see [ 135 ]) (see Algorithm 3 ((Note. The tabu search procedure is in the role of the self-contained algorithm.))), which, in turn, uses a self-contained tabu search algorithm—the “kernel” tabu search procedure. It is this procedure that directly improves a given solution. This procedure is thus in the role of the search intensification mechanism, while the mutation procedure is responsible for the diversification of the search process. It can be seen that the structure of the individual hierarchical levels of the HITS algorithm is quite simple, but the overall efficacy of the resulting multi-level algorithm increases significantly, which is demonstrated by the computational experiments. Of course, the run time increases as well, but this is compensated by the higher quality of the final results.

The interesting analogy between intensification and diversification (on the one side) and the phenomenon of entropy (on the other side) can be perceived. Indeed, the intensification process can be thought of as a process of the decrease of the entropy; on the other hand, diversification could be viewed as the increase of the entropy. Actually, the similar processes are seen in the open real physical systems. An example is the process of evolution of stars, where formation (birth) of the stars (along with the planets, organic matter, etc.) can be linked to the apparent decrease of the entropy, while the death of the stars (supernovae) may be associated with the significant increase of the entropy.

The self-contained tabu search procedure (for a more detailed description of the principles of TS algorithms, the reader is referred to [ 136 ]) iteratively analyses the neighbourhood of the current solution p (in our case— Θ 2 ( p ) ) and performs the non-prohibited move that most improves the value of the objective function. If there are no improving moves, then the one that least degrades the value of the objective function is accepted. In order to eliminate search cycles, the return to recently visited solutions is disabled for a specified period. The list of prohibitions—the tabu list T —is implemented as a two-dimensional matrix of size n × n . In this case, the entry t i j stores the sum of the number of the current iteration and the tabu tenure, h ; in this way, this value indicates from which iteration the i th and j th elements of a given solution can be again interchanged. The value of the parameter h depends on the problem size, n , and is chosen to be equal to 0.3 n . The tabu status is ignored at random moments with a very low probability α ( α ≤ 0.05 ). This allows to slightly increase the number of non-tabu solutions and not to limit the available search directions too much. The tabu condition is also ignored when the aspiration criterion is met, i.e., the current obtained solution is better than the best so far found solution. The resulting formal tabu and aspiration criteria are thus as follows:

t a b u _ c r i t e r i o n ( i , j ) = { TRUE , ( t i j ≥ q )   and   ( ς ≥ α )   and   ( H T ( ( z ( p ) + Δ ( p i , j ,   p ) ) mod   H a s h S i z e ) = TRUE ) FALSE , otherwise , a s p i r a t i o n _ c r i t e r i o n ( i , j ) = { TRUE , z ( p ) + Δ ( p i , j ,   p ) < z • FALSE , otherwise , where i , j are the current elements’ indices, q denotes the current iteration number, ς is a (pseudo-)random real number within the interval [ 0 ,   1 ] , and z • denotes the best so far found value of the objective function. H T denotes the hash table, which is simply a one-dimensional array, and H a s h S i z e is the capacity of the hash table.

In addition, our tabu search procedure uses a so-called secondary memory Γ [ 137 ] to avoid stagnation manifestations during the search process. The purpose of this memory is to gather high-quality solutions, which although are rated as very good, but are not chosen. In particular, the secondary memory includes solutions left “second” after the exploration of the neighbourhood Θ 2 . So, if the best found solution does not change more than ⌊ β τ ⌋ iterations, then the tabu list is cleared and the search is restarted from one of the “second” solutions in the secondary memory Γ (here, τ denotes the number of iterations of the TS procedure, and β is a so-called idle iterations limit factor such that 1 ≤ ⌊ β τ ⌋ ≤ τ ). The TS procedure is completed as soon as the total number of iterations, τ , has been performed.

The time complexity of the TS algorithm is proportional to O ( n 2 ) for the reason that the exploration of the neighbourhood Θ 2 requires n ( n − 1 ) 2 operations and also one needs to recalculate the differences of the objective function after each transition from a given solution to the new one.

The pseudo-code of the tabu search algorithm is shown in Algorithm 4 (Notes. (1) The immediate if function iif( x , y 1 , y 2 ) returns y 1 if x = TRUE , otherwise it returns y 2 . (2) The function random() returns a pseudo-random number uniformly distributed in [ 0 ,   1 ] . (3) The function random( x 1 , x 2 ) returns a pseudo-random number in [ x 1 ,   x 2 ] . (4) The values of the matrix Ξ are recalculated according to the Formula (9). (5) β denotes a random access parameter (we used β = 0.8 ).).

3.4.2. Mutation

Each solution found by the tabu search algorithm is subject to perturbation in the mutation procedure. Remind that formally the mutation process can be defined by the use of the operator φ : Π n → Π n . Thus, if p ~ = φ ( p ) , then p ~ ∈ Π n , p ~ ≠ p . More formalized operator can be described as follows: φ ( p ,   ξ ) : Π n × N → Π n , which transforms the current solution p to the new solution p ~ such that δ ( p ,   p ~ ) = δ ( p ,   φ ( p ) ) = ξ . In this way, 100 ξ n per cent elements of the solution are affected. The parameter ξ ( 2 ≤ ξ ≤ n ) regulates the strength of mutation and is referred to as a mutation rate. (In our algorithm, ξ = ⌊ 0.2 n ⌋ .) It is clear that for any p , p ~ (such that p ≠ p ~ , δ ( p ,   p ~ ) = ξ ) there always exists a sequence of distinct integers u 1 , u 2 , … ,   u ξ such that p ~ = ( ( ( p u 1 ,   u 2 ) u 2 ,   u 3 ) … ) u ξ − 1 ,   u ξ .

Choosing the right value of the mutation rate, ξ , plays a very important role in the mutation procedure and the HITS algorithm and, at the same time, the whole genetic algorithm. A proper compromise between two extreme cases must be achieved: (1) the value of ξ is (very) low (close to 0 ); (2) the value of ξ is (very) high (close to n ). In the first case, the mutation would not guarantee that the obtained mutated solution is “far” away enough from a given solution, which would lead to cycling search trajectories. In the second case, useful accumulated information would be lost and the algorithm would become very similar to a crude random multi-start method.

It should be stressed that the mutation processes are quite different from the crossover procedures. Mutation processes are by their nature purely random processes. Whilst crossover procedures only recombine the genetic code contained in the parents, the mutation processes generate, in essence, new information that hadn’t existed in predecessors earlier. It is the mutation process that implicitly is a true creative process and potentially produces a real novelty. In our work, twelve different mutation procedures and their modifications were proposed and tested.

3.4.2.1. Mutation Based on Random Pairwise Interchanges (M1)

In the beginning, the sequence r = ( r ( 1 ) ,   r ( 2 ) ,   … ,   r ( ξ ) ) of random integers r ( i ) ∈ { 1 ,   … . ,   n } is generated. Then, we start with the pair ( r ( 1 ) ,   r ( 2 ) ) , and the elements p ( r ( 1 ) ) , p ( r ( 2 ) ) are interchanged. Then, we exchange the elements p ( r ( 2 ) ) , p ( r ( 3 ) ) , and so on. This is repeated ξ − 1 times, where ξ is the value of the mutation rate defined by the algorithm’s user. The result of the mutation procedure is thus the solution p ~ satisfying the conditions: p ~ ∈ Π n , δ ( p ,   p ~ ) = ξ (see Figure 3 ).

Illustration of the mutation procedure ( n = 9 , ξ = 5 ) (The mutation process steps are as follows: (1) element 3 is interchanged with element 4; (2) element 3 (in position 3) is interchanged with element 1; (3) element 3 (in position 4) is interchanged with element 8; (4) element 3 (in position 6) is interchanged with element 5 (element 3 is eventually in position 8)).

On the basis of the random pairwise interchanges, other modified mutation procedures can be developed [ 138 ].

3.4.2.2. Random Pairwise Interchanges Using Random Key (M2)

In this case, the mutation process consists of two basic steps: (1) random pairwise interchanges; (2) shuffling the interchanged elements using a random key. A random key, r k , is a list of random indices of size ξ — r k ( 1 ) ,   r k ( 2 ) ,   … ,   r k ( ξ ) . The random key values simply determine which elements are again interchanged. The intention is to get a more “deeply” mutated solution and avoid returning to previously visited solutions.

3.4.2.3. Mutation Using Opposite Values (M3)

In this mutation procedure, the position’s index, let’s say k , is randomly determined. Then, the element e = p ( k ) is replaced by the following opposite value: o = ( ( p ( k ) + n 2 − 1 )   mod   n ) + 1 , where mod denotes the modulo operation. After this replacement, the solution element that was previously equal to o must also be changed. After both changes, p ( k ) becomes equal to o , p ( l ) —equal to e ; l indicates the element which was equal to o . The process is repeated ξ 2 times, where ξ is the muation rate.

3.4.2.4. Distance-Based Mutation (M4)

In this procedure, the indices of the pairs of elements to be interchanged are generated in such a way that the “distance” ( d ) between those indices is as large as possible. The following formula for generating the indices k 1 ,   k 2 ,   … ,   k ξ is used: k l = ⌊ ( ( d q l + ς − 1 )   mod   n ) + 1 ⌋ , here d = n ξ , ς —(pseudo)random real number from the interval [ 0 ,   1 ] , q l = ( q l − 1   mod   n ) + 1 , l = 1 , 2 , … ,   ξ ; the initial value q 0 is a random integer from the interval [ 1 , n ] .

3.4.2.5. Modified Random Pairwise Interchanges—Variant I (M5)

This is similar to the random pairwise interchanges. The sequence of random real-coded values from the interval [ 0 ,     1 ] is generated. The generated numbers along with their corresponding indices—known as smallest positive values—are sorted in the ascending order. These values, in particular, determine the elements to be interchanged.

3.4.2.6. Modified Random Pairwise Interchanges—Variant II (M6)

The list of random indices is obtained by directly generating random integers from the interval [ 1 ,   n ] . The integers may duplicate each other. To avoid duplications, the integers are sorted according to the ascending order. Indices corresponding to the sorted numbers indicate the elements that are to be interchanged.

3.4.2.7. Combined Mutation (M7)

This mutation procedure consists of two combined mutation procedures. Initially, the list of indices of the pairs of elements to be interchanged is constructed (see Section 3.4.2.6 ). The selected elements are then changed using opposite values (see Section 3.4.2.3 ).

3.4.2.8. Greedy Adaptive Search Based Mutation (M8)

The basic principle of this mutation procedure is that the solution is disintegrated in some way, and then reconstructed. The mutation process consists of two steps: (1) disintegration of the solution, which is random; (2) reconstruction of the solution, which is greedy-deterministic. In the first step, ξ elements are disregarded. In the second step, a greedy constructive algorithm is applied, which tries to find the best possible solution out of ξ ! available options. The value of ξ in this case should be quite small to prevent large increase in the run time of the mutation procedure. This mutation procedure (and also other procedures described below) are no longer problem-independent as the problem domain-specific knowledge is taken into account.

3.4.2.9. Greedy Randomized Adaptive Search Based Mutation (M9)

This mutation procedure resembles the one described above. The difference is that a greedy randomized adaptive search procedure (GRASP) [ 70 ] is used in the partial solution reconstruction phase to obtain an improved solution.

3.4.2.10. Randomized Local Search Based Mutation—Variant I (M10)

In this case, quick procedure based on random pairwise interchanges is initially performed (see Section 3.4.2.1 ). Then, a set of randomly selected elements is formed. A local search is then performed using the constructed set, i.e., only transitions between solutions that improve the value of the objective function are accepted.

3.4.2.11. Randomized Local Search Based Mutation—Variant II (M11)

This mutation variant is similar to the previous randomized local search variant, except that the randomly constructed neighbourhood is fully explored in a systematic way. Again, only improving transitions between solutions are accepted.

3.4.2.12. Randomized Tabu Search Based Mutation (M12)

Let p ( 1 ) = argmin i = 1 ,   … ,   n − 1 ,   j = i + 1 ,   … ,   n ,   m o v e _ a c c e p t a n c e _ c r i t e r i o n ( i , j ) = TRUE { z ( p i , j ) } and

p ( 2 ) = argmin i = 1 ,   … ,   n − 1 ,   j = i + 1 ,   … ,   n ,   m o v e _ a c c e p t a n c e _ c r i t e r i o n ( i , j ) = TRUE { z ( p i , j ) | p i , j ≠ p ( 1 ) } , where: m o v e _ a c c e p t a n c e _ c r i t e r i o n ( i , j ) = { TRUE , ( t a b u _ c r i t e r i o n ( i , j ) = FALSE )   or   ( a s p i r a t i o n _ c r i t e r i o n ( i , j ) = TRUE ) FALSE , otherwise , t a b u _ c r i t e r i o n ( i , j ) = { TRUE , ( t i j ≥ q )   and   ( ς ≥ α ) FALSE , otherwise ,

a s p i r a t i o n _ c r i t e r i o n ( i , j ) = { TRUE , z ( p ) + Δ ( p i , j ,   p ) < z • FALSE , otherwise , q denotes the current iteration number, ς is a (pseudo-)random number within the interval [ 0 ,   1 ] , α denotes the randomization coefficient, z • denotes the best so far value of the objective function. Then, in the randomized tabu search procedure, the best achieved solution (“winner solution”) p ( 1 ) is accepted with the probability γ , meanwhile the second solution p ( 2 ) is chosen with the probability 1 − γ (note that, in the case of γ = 1 , we get the standard (deterministic) tabu search.) In our algorithm, we use γ = 0.2 . So, the central idea of the randomized tabu search is just this quasi-random mixing between the “winner solution” and “next to the winner solution” in the course of the tabu search process. Based on the extensive experimentation, we found out that this type of mutation is the most promising mutation procedure among all the procedures examined. The explanation would be that this type of mutation rather is more gentle, moderate and controlled than the other mutation procedures.

In the end, note that the computational complexity of all our mutation procedures is proportional to O ( ξ n 2 ) . This is due to the fact that our mutation procedures recalculate the differences of the objective function (i.e., the values of the matrix Ξ ) approximately ξ times (see Algorithm 5 (Note. The values of the matrix Ξ are recalculated using Equation (9).)). So, the smaller the value of ξ , the faster is the mutation procedure. Also, note that the difference matrix Ξ is (permanently) stored in a RAM (operating memory), so there is no need to calculate the differences of the objective function from scratch.

3.4.3. Candidate Acceptance

Regarding the candidate solution acceptance rule, we always choose the most recently (newly) found improved solution (the latest result of the HITS (or TS) algorithm) instead of the overall best found solution. Such an approach is thought to allow to explore potentially larger regions of the solution space.

3.5. Population Replacement

For the population replacement, a modified rule is used to respect not only the quality of the solutions, but also the difference (distance) between solutions.

We have, in particular, implemented an extended variant of the well-known “ μ + λ “ update rule [ 139 ]. The new advanced replacement rule is denoted as “ μ + λ , ε ”. (This rule is also used for the initial population construction (see Section 3.1 ).) We remind that if the minimum mutual distance between population members and the new obtained offspring is less than the distance threshold, D T , then the offspring is omitted. The only exception is the case where the offspring appears better than the best population member. Otherwise, the offspring enters the current population, but only under condition that it is better than the worst population member. The worst individual is removed in this case. This our replacement strategy ensures that only individuals that are diverse enough survive for the further generations.

There are a few replacement variations (options), depending on the parameter R e p V a r . If R e p V a r = 1 , then exactly the above replacement strategy is adopted. If R e p V a r = 2 , then the new offspring replaces the best member of the current population if the offspring is better than the best population individual. If the offspring is worse than the best individual, then R e p V a r = 2 is identical to R e p V a r = 3 . If R e p V a r = 3 , then the offspring replaces the worst member of the population, ignoring the fitness of the worst individual. The minimum distance criterion must be taken into account.

3.6. Population Restart

The important feature of our genetic algorithm is the use of the population restart mechanism to try to avoid the premature convergence and search stagnation. The restart process is triggered in the situations where the solutions of the population do not change at all for some number of consecutive generations. This can be operationalized by the use of a priori parameter called an idle generations limit, L i d l e _ g e n   , where L i d l e _ g e n = max { L m i n ,   ⌊ ω N g e n ⌋ } , here L m i n is a constant (we use L m i n = 3 ), ω is to denote a stagnation control parameter and N g e n is the total number of generations of the genetic algorithm. (The standard value of ω is equal to 0.05 .) The restart itself is performed by applying a so-called multi-mutation, where the mutation process is applied to all the members of the stagnated population. Such approach is preferred to the complete destroying of the population, which seems to be too aggressive.

4. Computational Experiments

Our genetic-hierarchical algorithm is implemented by using C# programming language. The computational experiments have been carried out on a 3.1 GHz personal computer running Windows 7 Enterprise. The CPU model is an Intel Core i5-3450.

The algorithm is tested on the small-, medium-and large-scaled QAP benchmark instances of sizes from n = 10 to n = 128 . Most instances are from the online QAP library QAPLIB [ 29 ]. Other instances are from [ 14 , 19 ] (see also http:/mistic.heig-vd.ch/taillard/problemes.dir/qap.dir/qap.html).

In particular, the following benchmark instances taken from QAPLIB were examined:

• (a) random, unstructured instances (these instances are denoted by: rou20, tai10a, tai12a, tai15a, tai17a, tai20a, tai25a, tai30a, tai35a, tai40a, tai50a, tai60a, tai80a, tai100a);
• (b) randomly generated, grid-based instances (they are denoted by: had20, nug30, scr20, sko42, sko49, sko56, sko64, sko72, sko81, sko90, sko100a..sko100f, tho30, tho40, wil50, wil100);
• (c) real-life, structured instances from practical applications (denoted by: chr25a, els19, esc32a..esc32h, esc64a, esc128, kra30a, kra30b, ste36a. ste36c, tai64c);
• (d) real-life like (pseudo-random), structured instances (denoted by: tai10b, tai12b, tai15b, tai20b, tai25b, tai30b, tai35b, tai40b, tai50b, tai60b, tai80b, tai100b).
• (e) instances with known optimal solutions (denoted by: lipa20a, lipa20b, lipa30a, lipa30b, lipa40a, lipa40b, lipa50a, lipa50b, lipa60a, lipa60b, lipa70a, lipa70b, lipa80a, lipa80b, lipa90a, lipa90b).

In addition, the instances introduced by de Carvalho and Rahmann [ 19 ] are investigated. These instances are extremely difficult to solve. They are denoted by bl36, bl49, bl64, bl81, bl100 (aka. border length minimization instances) and ci36, ci49, ci64, ci81, ci100 (aka. conflict index minimization instances).

We also tested the instances dre15, dre18, dre21, dre24, dre28, dre30, dre42, dre56, dre72, dre90, tai27e1, tai27e2, tai27e3, tai45e1, tai45e2, tai45e3, tai75e1, tai75e2, tai75e3 proposed by Drezner and Taillard in [ 14 ].

In the initial computational experiments, we used the following “learning set” of the QAP benchmark instances of sizes from n = 35 to n = 70 : bl49, bl64, ci49, ci64, dre42, dre56, lipa70a, lipa70b, sko56, sko64, tai35a, tai35b, tai40a, tai40b, tai45e1, tai50a, tai50b, tai60a, tai60b, wil50. These particular instances were chosen based on our preliminary experience.

As a performance criterion, we adopt the average relative percentage deviation ( θ ¯ ) of the yielded solutions from the best known solution (BKS). It is calculated by the following formula: θ ¯ = 100 ( z ¯ − z b k v ) z b k v [ % ] , where z ¯ is the average objective function value over 10 runs of the algorithm, while z b k v denotes the best known value (BKV) of the objective function that corresponds to the BKS (BKVs are from [ 14 , 29 , 86 ]).

At every run, the algorithm is applied to the particular instance. Each time, the algorithm starts from a new random initial population. The algorithm is terminated if either the maximum number of generations, N g e n , has been reached or the best known solution has been achieved.

In the experiments, the goal was to empirically test the performance of the basic setup of our algorithm and also its various variants in terms of the quality of solutions and the run time of the algorithm. To do so, we have identified some essential algorithm’s components (ingredients) (namely, “initial population”, “selection”, “crossover”, “local improvement (hierarchical tabu search)”, “mutation”, “population replacement”) to reveal their influence on the efficiency of GHA and to “synthesize” the preferable fine-tuned architecture of the algorithm. The following combination of the particular options (parameters) related to these components is declared as the basic version of GHA: { I n i t P o p V a r = 1 , P S = 10 , N g e n = 100 , σ = 1.5 , C r o s s V a r = “ 1 PX ” , τ = 20 , Q h i e r = 2 8 = 256 , M u t V a r = “ M 1 ” , R e p V a r = 1 } ; here, Q h i e r denotes the total cumulative number of hierarchical iterations ( Q h i e r = Q ( 0 ) × Q ( 1 ) × Q ( 2 ) × Q ( 3 ) × Q ( 4 ) × Q ( 5 ) × Q ( 6 ) × Q ( 7 ) ), Q ( 0 ) ,   … ,   Q ( 7 ) denote, respectively, the corresponding numbers of iterations of the 0th-level, …, 7th-level hierarchical iterated tabu search algorithms. The prescribed default values of the control parameters corresponding to the basic version of GHA are shown in Table 1 and Table 2 . (These values were later over-ridden in particular separate experiments).

Values of the control parameters of the basic version of GHA used in the experiments.

Standard values of the control parameters of the hierarchical iterated tabu search algorithm.

† Q ( 0 ) = Q ( 1 ) = Q ( 2 ) = Q ( 3 ) = Q ( 4 ) = Q ( 5 ) = Q ( 6 ) = Q ( 7 ) = 2 .

In the initial experiments, twelve crossover procedures (1PX, 2PX, UX, SX, PMX, SPX1, SPX2, SPX3, CX, COHX, MPX, UNIVX) have been compared against each other. The obtained results (presented in Table 3 ) demonstrate that our proposed universal crossover (UNIVX) with the tuned control parameters yields the most promising results.

Comparison of crossover procedures.

Notes. Time denotes the average CPU time per one run. In all cases, the first mutation variant (M1) is used.

In the further experiments, the different mutation procedures (M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M11, M12) were examined. This time, we have found out that the randomized tabu search based mutation is clearly the best among the all tested mutation variants (see Table 4 ).

Comparison of mutation procedures.

Note. In all cases, the 1PX crossover is used.

Further, we were interested in how various options (configurations) of the initial population construction affect the performance of the genetic-hierarchical algorithm. The particular separate configurations differ with respect to the option of the population construction ( I n i t P o p V a r ), the size of pre-initial population ( P S ′ ), as well as the number of TS iterations during the population initialization ( τ ′ ). In particular, the following variants were investigated: (1) I n i t P o p V a r = 1 , P S ′ = 10 , τ ′ = 20 ; (2) I n i t P o p V a r = 1 , P S ′ = 20 , τ ′ = 40 ; (3) I n i t P o p V a r = 1 , P S ′ = 40 , τ ′ = 80 ; (4) I n i t P o p V a r = 1 , P S ′ = 100 , τ ′ = 200 ; (5) I n i t P o p V a r = 2 , P S ′ = 10 , τ ′ = 20 ; (6) I n i t P o p V a r = 2 , P S ′ = 20 , τ ′ = 40 ; (7) I n i t P o p V a r = 2 , P S ′ = 40 , τ ′ = 80 ; (8) I n i t P o p V a r = 2 , P S ′ = 100 , τ ′ = 200 ; (9) I n i t P o p V a r = 3 , P S ′ = 10 , τ ′ = 20 ; (10) I n i t P o p V a r = 3 , P S ′ = 20 , τ ′ = 40 ; (11) I n i t P o p V a r = 3 , P S ′ = 40 , τ ′ = 80 ; (12) I n i t P o p V a r = 3 , P S ′ = 100 , τ ′ = 200 .

We have observed that maintaining the higher quality initial populations, in general, allows to significantly increase the overall efficiency of GHA when comparing to the lower quality initial populations (see Table 5 ).

Comparison of different variants of the initial population construction.

Note. In all cases, the UNIVX crossover and the twelfth mutation variant (M12) are used.

Additionally, we experimented with some few population replacement options. The particular population replacement variants are as follows: (1) P S ′ = 10 , τ ′ = 20 , R e p V a r = 1 ; (2) P S ′ = 10 , τ ′ = 20 , R e p V a r = 2 ; (3) P S ′ = 10 , τ ′ = 20 , R e p V a r = 3 ; (4) P S ′ = 20 , τ ′ = 40 , R e p V a r = 1 ; (5) P S ′ = 20 , τ ′ = 40 , R e p V a r = 2 ; (6) P S ′ = 20 , τ ′ = 40 , R e p V a r = 3 ; (7) P S ′ = 40 , τ ′ = 80 , R e p V a r = 1 ; (8) P S ′ = 40 , τ ′ = 80 , R e p V a r = 2 ; (9) P S ′ = 40 , τ ′ = 80 , R e p V a r = 3 ; (10) P S ′ = 100 , τ ′ = 200 , R e p V a r = 1 ; (11) P S ′ = 100 , τ ′ = 200 , R e p V a r = 2 ; (12) P S ′ = 100 , τ ′ = 200 , R e p V a r = 3 .

It was observed that the aggressive strategy of replacement of the best population member ( R e p V a r = 2 ) seems to be superior to other options (see Table 6 ). Further, more extensive experiments are required to strengthen this conjecture.

Comparison of different variants of population replacement.

Note. In all cases, the UNIVX crossover and the mutation variant M12 are used. Also, the initial population construction option I n i t P o p V a r = 1 is used.

In addition, we have tested some other different scenarios (regimes) in order to unveil some possible tendencies of the behaviour of the HITS algorithm. The following scenarios were investigated: (1) scenario of “quick search”: small value of τ —small value of Q h i e r ; (2) scenario of “diversified quick search”: small value of τ —large value of Q h i e r ; (3) scenario of “extensive search”: large value of τ —small value of Q h i e r ; (4) scenario of “diversified extensive search”: large value of τ —large value of Q h i e r . Note that, in these scenarios, the number of generations of GHA should be accordingly balanced in order to stay within the fixed run time. The corresponding values of the control parameters are as follows: (1-a) τ = 20 , Q h i e r = 2 8 = 256 , N g e n = 100 ; (1-b) τ = 25 , Q h i e r = 2 8 = 256 , N g e n = 80 ; (1-c) τ = 40 , Q h i e r = 2 8 = 256 , N g e n = 50 ; (2-a) τ = 50 , Q h i e r = 2 8 = 256 , N g e n = 40 ; (2-b) τ = 10 , Q h i e r = 2 9 = 512 , N g e n = 100 ; (2-c) τ = 20 , Q h i e r = 2 9 = 512 , N g e n = 50 ; (3-a) τ = 10 , Q h i e r = 5 × 2 7 = 640 , N g e n = 80 ; (3-b) τ = 20 , Q h i e r = 5 × 2 7 = 640 , N g e n = 40 ; (3-c) τ = 10 , Q h i e r = 2 10 = 1024 , N g e n = 50 ; (4-a) τ = 20 , Q h i e r = 2 10 = 1024 , N g e n = 25 ; (4-b) τ = 10 , Q h i e r = 10 × 2 7 = 1280 , N g e n = 40 ; (4-c) τ = 20 , Q h i e r = 10 × 2 7 = 1280 , N g e n = 20 .

The results of the experiments (see Table 7 ) demonstrate that the scenario of diversified extensive search is obviously preferable to other examined scenarios.

Comparison of different variants of the hierarchical iterated tabu search algorithm.

Note. In all cases, the UNIVX crossover and the mutation variant M12 are used. Also, the initial population construction option I n i t P o p V a r = 1 and population replacement option R e p V a r = 2 are used.

Additional scenarios have been examined to reveal the reaction of GHA when extensively increasing the cumulative number of iterations of the hierarchical iterated tabu search algorithm— Q h i e r . The computational budget is not constant (“balanced”) anymore, but grows as the value of Q h i e r increases. The following scenarios have been tried: (1) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 5 × 2 7 = 640 ; (2) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 6 × 2 7 = 768 ; (3) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 7 × 2 7 = 896 ; (4) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 8 × 2 7 = 1024 ; (5) P S = 20 , P S ′ = 300 , τ ′ = 300 , Q h i e r = 5 × 2 7 = 640 ; (6) P S = 20 , P S ′ = 300 , τ ′ = 300 , Q h i e r = 6 × 2 7 = 768 ; (7) P S = 20 , P S ′ = 300 , τ ′ = 300 , Q h i e r = 7 × 2 7 = 896 ; (8) P S = 20 , P S ′ = 300 , τ ′ = 300 , Q h i e r = 8 × 2 7 = 1024 ; (9) P S = 10 , P S ′ = 200 , τ ′ = 400 , Q h i e r = 7 × 2 7 = 896 ; (10) P S = 10 , P S ′ = 200 , τ ′ = 400 , Q h i e r = 7 × 2 7 = 1024 ; (11) P S = 20 , P S ′ = 400 , τ ′ = 400 , Q h i e r = 7 × 2 7 = 896 ; (12) P S = 20 , P S ′ = 400 , τ ′ = 400 , Q h i e r = 7 × 2 7 = 1024 ; here, τ ′ denotes the number of TS iterations during the construction of the initial population.

The results confirm that, as expected, there exists a clear correlation between the number of improving iterations (the number of TS/HITS iterations) and the quality of the obtained solutions (see Table 8 ).

Results of the experiments with the increased number of iterations of the hierarchical iterated tabu search algorithm.

Notes. In all cases, the UNIVX crossover and the mutation variant M12 are used. Also, the options I n i t P o p V a r = 3 , R e p V a r = 2 are used.

To have a reflection of the obtained results from a different perspective—in particular, a demonstration of the stability and robustness properties of our algorithm—we have constructed histograms of the frequency of the objective function values for one of the most difficult instances of the “learning set”—bl64 (see Figure A1 in the “ Appendix A ” Section). In fact, we have created the histograms of the frequency of the average percentage deviation, θ ¯ , over 10 algorithm runs within the interval [ 0.0 ,   1.0 ) , where 0.0 stands for zero deviation and 1.0 means the maximum possible deviation. (Note that the average deviation never exceeded 1.0 for the instance bl64 (see Table 8 ).)

(Regarding the selection factor, σ , the obtained results are quite “flat” and not statistically significant, so they are omitted).

On the whole, we have found the best known solutions in the 9191 cases (runs) out of 14400 cases ( 64 % of cases). The BKS was found at least once for all examined instances. The cumulative average percentage deviation is equal to 1.452 % and the cumulative average CPU time per run is equal to approximately 65.9 s. The average deviation is less than 0.5 in 73 % of cases, while the average deviation is less than 1.0 in 89 % of cases. 14 instances (ci49, ci64, dre42, lipa70a, lipa70b, sko56, sko64, tai35a, tai35b, tai40b, tai45e1, tai50b, tai60b, wil50) were solved to pseudo-optimality in more than 300 runs.

After experimenting with the “learning set” of instances, the other instances (the “testing set” of instances) were examined using the fine-tuned parameters in order to find out how quickly the genetic-hierarchical algorithm converges to the best known/optimal solutions. The obtained results are presented in Table 9 . It can be seen that all tested instances ( 88 instances) are solved to pseudo-optimality within extremely small computation time.

Results of GHA for the set of 88 instances of QAPLIB [ 14 , 19 , 29 ].

Note. Time denotes the average CPU time per one run.

We have also compared our algorithm with the memetic algorithm (MA) proposed be Benlic and Hao [ 78 ], which is most likely the best so far heuristic algorithm for the QAP, to the best of our knowledge. The results of comparison of the algorithms are presented in Table 10 , Table 11 and Table 12 . It seems that our genetic-hierarchical algorithm outperforms MA. Additionally, we used the genetic algorithms by Drezner et al. [ 14 ] and Drezner and Misevičius [ 86 ] in the further comparison (see Table 13 , Table 14 , Table 15 and Table 16 ). Again, our algorithm compares favourably to both the algorithm by Drezner et al. as well as Drezner and Misevičius.

Comparative results between GHA and memetic algorithm (MA) [ 78 ] (part I).

Comparative results between GHA and memetic algorithm (MA) [ 78 ] (part II).

Notes. Time denotes the average CPU time per one run. In parentheses, we present the number of times that the BKS has been found. The best known value for tai100a is from [ 140 ].

Comparative results between GHA and memetic algorithm (MA) [ 78 ] (part III).

Comparative results between GHA and hybrid genetic algorithm (HGA) [ 14 ] (part I).

Note. Time denotes the average CPU time per one run. In parentheses, we present the number of times that the BKS has been found.

Comparative results between GHA and hybrid genetic algorithm (HGA) [ 14 ] (part II).

Comparative results between GHA and hybrid genetic algorithm with differential improvement (HGA-DI) [ 86 ] (part I).

Comparative results between GHA and hybrid genetic algorithm with differential improvement (HGA-DI) [ 86 ] (part II).

5. Concluding Remarks

In this paper, we have presented the hybrid genetic-hierarchical algorithm for the solution of the quadratic assignment problem. The key feature of the proposed algorithm is that the genetic algorithm is hybridized with the hierarchicity-based (self-similar) iterated tabu search algorithm, which serves as a powerful local optimizer of the offspring solutions produced by the crossover operator.

The algorithm was examined on the QAP benchmark instances of various sizes and complexity. The results obtained from the experiments demonstrate the excellent performance of the genetic-hierarchical algorithm. Our algorithm seems to outperform other state-of-the-art heuristic algorithms for many examined QAP instances or is at least very much competitive with them. A more pronounced improvement in the quality of the results might be achieved by a thorough calibration of the algorithm’s parameters.

The following are some possible future research directions: balancing of the number of tabu search iterations and the number of hierarchical iterated tabu search iterations, as well as the number of hierarchical levels; extensive experimental analysis of the particular components and configurations of the genetic-hierarchical algorithm; designing and implementing a multi-level hierarchical (master-slave) genetic algorithm.

Histograms of the frequency of the objective function values for the instance bl64 for different examined scenarios: ( a ) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 640 ; ( b ) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 768 ; ( c ) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 896 ; ( d ) P S = 10 , P S ′ = 150 , τ ′ = 300 , Q h i e r = 1024 . The histograms are developed in such a way that the frequency of the average deviation is visualized over 10 discrete sub-intervals of the interval [ 0.0 ,     : [ 0.0 ,   0.1 ) ; [ 0.1 ,   0.2 ) ; [ 0.2 ,   0.3 ) ; …; [ 0.9 ,   1.0 ) . It can be seen that the average deviation from (pseudo-)optimal solutions stably decreases by increasing the number of search iterations ( Q h i e r ).

Author Contributions

The proposed algorithm was designed and implemented by A.M. All sections and experiments were described by both authors. All authors have read and agreed to the published version of the manuscript.

This research was funded by the Faculty of Informatics of Kaunas University of Technology.

Institutional Review Board Statement

Informed consent statement, data availability statement, conflicts of interest.

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.