Bird Mating Optimizer for Combinatorial Optimization Problems

The bird mating optimizer is a new metaheuristic algorithm that was originally proposed to solve continuous optimization problems with a very promising performance. However, the algorithm has not yet been applied for solving combinatorial optimization problems. Thus, the formulation may not be able to generate a discrete feasible solution. Many continuous algorithms used random-key representation to represent the discrete solution using real numbers or a discrete variant of the algorithm is used to deal with the discrete solution of the problem. However, there is no evidence which methodology is better for solving combinatorial optimization problems. Therefore, this work proposes two variants of bird mating optimizer (random-key bird mating optimizer and the discrete bird mating optimizer), to identify which one is more efficient in solving combinatorial optimization problems. In the first one, we use a random-key encoding scheme, whilst, in the later one, we use crossover (multi-parent) and mutation operators to combine the components of the selected parents to generate new broods. The performance of these algorithms is tested on the travelling salesman problem and berth allocation problem, and are compared with the results of two well-known optimization algorithms: Genetic Algorithm and Particle Swarm Optimization. Experimental results show that the discrete bird mating optimizer is more efficient than the others on all tested benchmark instances. Indeed, it is able to attain the best-known results in some of the BAP benchmark instances. These indicate the applicability and the effectiveness of the proposed discrete bird mating optimizer in solving combinatorial optimization problems.


I. INTRODUCTION
Combinatorial optimization problems arise in various aspects of computer science and other areas such as artificial intelligence, operational research and electronic commerce. The combinatorial optimization problem can be defined as the need to find an optimal arrangement, grouping, or ordering for a given set of discrete variables [1], [2]. Examples of such problems include scheduling public transportation, travelling salesman problems, university educational timetabling, and berth allocation problems [3]. Several optimization techniques have been proposed for solving combinatorial optimization problems efficiently [1], [4]- [7], [74], [75].
The associate editor coordinating the review of this manuscript and approving it for publication was Fuhui Zhou .
These techniques can be classified as either exact methods or approximate methods. Exact methods such as dynamic programming [8], [9], branch and cut [10] and branch and price [11] are guaranteed to find an optimal solution in a finite time and systematically search the solution space [12]. However, due to the complexity of many combinatorial problems and many are known to be non-deterministic polynomial-time (NP)-hard problems, the time needed to solve them grows exponentially as the problem size grows linearly [13]. Therefore, many researchers have focused on approximate methods to address combinatorial optimization problems [13]- [15].
The use of approximate algorithms does not guarantee that the optimal solution will be found, but empirically they have often been shown to find a nearly optimal solution VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ within a reasonable amount of time [13]. Approximate algorithms are classified into two families: heuristics and metaheuristics. Heuristics, such as construction algorithms, are problem dependent, that is, they are designed to be applicable to a particular problem domain [16], [17]. Heuristic algorithms can be costly to implement and may not be suitable for a large variety of optimization problems [2], [17]. In contrast, metaheuristics are more general methodologies that can work effectively across different optimization problems. A metaheuristic algorithm is a high-level heuristic methodology that provides a high-level control strategy in order to improve exploration of the search space [18]. Examples of these algorithms include simulated annealing [19], [20], great deluge [21], tabu search [22], genetic algorithm (GA) [23]- [26], harmony search [27], scatter search [28], ant colony optimization [29], [30] and particle swarm optimization [31], [32].
Recently, a nature-inspired metaheuristic population-based algorithm, called the bird mating optimizer (BMO) has been proposed [14]. The BMO is inspired by the mating strategies of bird species during mating season in which male birds tend to mate with female birds in order to breed a new brood. The strength of the BMO lies in its ability to provide a good balance between exploration and exploitation [14]. The BMO overcomes the drawback of losing population diversity and getting trapped in local optima that is inherent in other methodologies by employing five updating strategies to move through the search space [14]. Furthermore, the BMO is able to exploit the best regions by utilizing a local search at each iteration. Therefore, BMO has more capability than the GA to effectively explore and exploit the search space and find the global solution. The algorithm has demonstrated effective performance across range of optimization problems [33]- [40]. In addition, it showed a competitive result comparing to other evolutionary algorithm such as classical evolutionary programming, fast evolutionary programming, classical evolutionary strategies, fast evolutionary strategies, genetic algorithm, particle swarm optimization, and group search optimizer [14].
However, BMO was originally proposed to address continuous optimization problems and, therefore, the application of BMO to combinatorial optimization problems has not yet been investigated. The original BMO uses mathematical formulations that combine the information of the selected solutions to generate a new one. Thus, standard BMO equations may not be able to generate a discrete feasible solution because the positions are discrete values in combinatorial optimization problems. Many continuous algorithms in the literature that have been applied to solve combinatorial problems use one of the two following methods. The first method is to apply a continuous algorithm using random-key representation to represent the discrete solution using real numbers, which can deal with the continuous algorithms [41]- [44]. The second method is to construct a discrete variant of the algorithm to deal with the discrete nature of the problem [15], [43], [45], [46]. However, there is no investigation in the literature on which method is better for solving combinatorial problems. The question is how can BMO be applied to solve combinatorial optimization problems and which method is better? To answer these questions, this work proposes two variants of BMO to address combinatorial optimization problems: the random-key bird mating optimizer (RKBMO) and discrete bird mating optimizer (DBMO). We also investigate the effectiveness of each method to identify the most promising one. Two well-known combinatorial optimization problems are chosen to test the proposed methods: The Travelling Salesman problem (TSP) (Reinelt 1991 instances) [47] and the Berth Allocation Problem (BAP) (Cordeau et al. 2005 instances) [48] Note this research study is different from [49], in which only RKBMO was applied and only a set of BAP instances were used and they are different from those the set of BAP used in this work. In addition, no analysis study was provided. Here we utilise two problem domains (TSP and BAP) and conduct deeper analysis.
The remainder of the paper is organized as follows: the basic BMO is presented in Section II. The problem descriptions are presented in Section III. The proposed RKBMO is presented in Section IV, followed by the proposed DBMO in Section V. In Section VI, the application of the proposed DBMO on the BAP and TSP is presented. Section VII summarizes and discusses the computational results. Finally, a conclusion and details of future work are presented in Section VIII.

II. THE BIRD MATING OPTIMIZER
The BMO is a population-based stochastic search technique that was proposed by [14] to address continuous optimization problems. The behavior of this algorithm is based on simulating the mating strategies of bird species during mating season. In this algorithm, the mating process of birds involves the use of three main operators to produce a new generation: two-parent mating, multi-parent mating, and mutation. Two-parent mating is where the two parents mate with each other to breed only one new brood, whereas multi-parent mating is when a parent mates with at least two other parents to breed one new brood. Mutation, on the other hand, is a mechanism where the female parent produces a new brood without the help of a male by modifying its own genes [14].
In the BMO, the population is represented by the bird society, and each bird represents a feasible solution to the problem at hand. There are five mating types that can generate a new population: Thus, a bird society is divided into males and females. The females are those birds that have the most promising genes in the society. The females are divided into two groups: parthenogenetic and polyandrous, while the males are divided into three groups: monogamous, polygynous and promiscuous [14].
Generally, in the original BMO, the assumption that each bird in the society can change its own strategy to another when a new generation is updated [14]. At each generation, the birds that have the highest fitness values are designated as parthenogenetic and polyandrous, while those that have worse fitness values are designated as promiscuous. The remaining birds in the society are considered as monogamous or polygynous birds. A more detailed explanation of each group is given below.
Parthenogenesis is the mating type in which the female bird can produce a brood without mating with a male. In this method, each female tries to produce her brood by modifying and changing her genes with a predefined rate. Each female bird in the parthenogenetic group produces a brood as presented in equation (1) [14]: where x (i) denotes the i th bird, x brood is the resultant brood, n is the problem dimension (number of genes), mcf p is the mutation control factor of parthenogenesis, r 1 denotes a random number between 0 and 1, and µ is the step size. Monogamy is a two-parent mating type in which a male tends to mate with only one female. Every male evaluates the quality of the females by using a probabilistic approach to select one of them to be its mate. Female birds with good genes have a higher probability of being selected. Equation 2 illustrates the process of producing a new brood from two selected parents [13]: where w is a time-varying weight to adjust the selected female, r is a 1 × d vector in which every element is a distributed random number between 0 and 1 and this random vector influences the corresponding element of ( x i − x ), n is the problem dimension, mfc denotes the mutation control factor which is distributed between 0 and 1, and u and l are the upper bounds and lower bounds of the elements, respectively.
From the first part of Equation (2), each male bird seeks to produce a good-quality brood by finding a worthy female with which to mix his genes in order to produce new genes. Then, the male bird tries to improve his produced brood by mutating one of the brood's genes with the probability of 1 − mcf .
Polygamy is a multi-parent mating process in which a male bird tries to produce a brood by mating with two or more females. The benefit of this multi-mating process is to produce a brood with better genes. In nature, a polygynous bird mating with several females produces a number of broods, but in the BMO, only one brood results from this mating process, where the brood's genes are a combination of the male's and multiple females' genes. After selecting the females by using a selection mechanism, each male bird mates with his selected females. The resultant brood is produced by the following process [13]: where n i is the number of selected female birds and x i j denotes the jth selected bird.
Promiscuity is a two-parent mating approach between two birds where they only mate once. This type of mating does not lead to a long relationship between the birds. This mating process indicates a chaotic social structure in which the male bird will never see the brood or the nest, and generally will not see the female for further mating activity. In promiscuity, the birds use a chaotic sequence method during the generations. However, each promiscuous bird behaves in the same way as a monogamous bird. In other words, in the BMO, promiscuous birds produce a new brood according to Equation (2).
Polyandry is a multi-parent mating system where a female bird seeks to mate with more than two males. The female performs a selection mechanism to select the males. Then, each female bird mates with her selected male birds. In the BMO, polyandrous birds produce a new brood according to Equation (3).
The procedure of the BMO is as follows [14]: Step 1 (Parameter Initialization): initialize the following BMO parameters: 1) Society size (SS): this refers to the number of birds in the society or the number of solutions in the population of the BMO; 2) The percentage of each group in the society: that is, the percentage of monogamous, polygynous, promiscuous, polyandrous, and parthenogenetic birds; 3) The number of mates (nm) for the polygynous and polyandrous birds: this refers to the number of mates that will participate in the multi-parent mating system practiced by the polygynous and polyandrous birds; 4) Mutation control factor (mcf): this is the probability of mutating the generated brood after each mating process; 5) Maximum number of generations (NG): this refers to the termination condition of the BMO or the number of iterations.
Step 2 (Society Initialization): randomly initialize a set of feasible solutions and add them to the society. Each solution is considered as a bird and is specified by a vector, with the length of n.
Step 3 (Society Evaluation): calculate the quality of each bird using an evaluation function.
Step 4 (Ranking): rank the birds in the society based on their quality.
Step 5 (Classification): divide the society into five groups of birds: parthenogenetic, polyandrous, monogamous, polygynous and promiscuous. The parthenogenetic and polyandrous groups are considered as female and the remaining groups are considered as male. Birds with top n1 qualities are designated as parthenogenetic and the next n2, n3, n4 and n5 are designated as polyandrous, monogamous, polygynous and promiscuous, respectively. The percentage of each group is defined in the parameter initialization step, as recommended by [14].
Step 6 (Breeding): each bird produces a new brood using its own pattern.
Step 7 (Replacement): if the quality of the brood is better than the quality of the birds in the society, then the brood replaces the bird, otherwise, the bird remains in the society and the brood is abandoned.
Step 8 (Termination Condition): repeat steps 4 to 7 until a predetermined number of generations (NG) is performed.
Step 9 (Report the Best): select the best-quality bird in the society as the best solution.

III. PROBLEM DESCRIPTION
In this study, the performance of the proposed DBMO is evaluated over two benchmark problems: The Travelling Salesman Problem (TSP) and Berth Allocation Problem (BAP). The following subsections describe these problems.

A. TRAVELLING SALESMAN PROBLEM
The TSP is a well-known combinatorial optimization problem that classified as a NP-hard problem, which means that it may take an infeasible computational time to solve it [50]- [52]. The TSP can be described as a search for the shortest path that visits each city once and only once, and finally return to the start city [53].

B. BERTH ALLOCATION PROBLEM
The BAP is a highly constrained combinatorial optimization problem which is hard to solve to optimality [48], [54]. The BAP can be classified into two types based on berth type and vessel arrival time. The berth type is considered as discrete if the quay is divided into a set of sections (berths) or continuous if not partitioned. The vessel arrival time is considered as dynamic if vessels can arrive at any time during the container operations with planning arrival time or static if all the vessels have to arrive in the harbor before the berth planning step begins. In this study, we focus on the discrete and dynamic version of the BAP [43]. This BAP deals with allocating vessels to berths in the harbor at the planned arrival time. More formally, the goal is to assign a berth for each vessel and a service time at the selected berth. The following assumptions are considered in the BAP [46]: 1) Each berth can serve only one vessel at a time; 2) Any vessel can be assigned to any berth with a given handling time taking into account that the handling time of a vessel can differ from one berth to another; 3) All vessels arrive at their berths before or after the berths' opening hours with a known arrival time; 4) When a vessel is moored in a berth, it remains there until all servicing activities have been completed. The objective of the BAP is to minimize the overall waiting time of all the vessels that need to be serviced in the harbor, which is calculated as an objective function as follows [48]: Random-key representation is a common technique that transforms a position in a continuous space and converts it onto a combinatorial space representation [41]. It uses a vector of real numbers to represent a solution in which each number is randomly generated in uniform [0, 1]. The combinatorial vector is composed of integers ordered according to the sequence of the real numbers in the first vector. This scheme is proposed by [41] and used in the literature to address different combinatorial optimization problems [44], [43], [55]. An example of this representation is shown in Fig. 1.  A. RKBMO FOR TSP Fig. 2 illustrates the steps of encoding and decoding process of generating a TSP solution using random keys. To encode the solution for BMO, we generate m random numbers from [0, 1]m for each bird in the society, where m is the size of the TSP instance. Each gene in the bird has an integer index (see Fig. 4). In the decoding process, we follow the steps as presented in [44], [41]: 1) Sort the random keys in ascending order. 2) Set the indexes of the sorted numbers as the city index array, i.e. the integer indexes correspond to the city index. The sequence of the generated real numbers represents the order of the genes in a bird. Thus, a bird contains a set of genes (i.e. city index).
During the search process, some numbers might be outside the feasible solution space. We handle this limitation by doing the following: where the solution is infeasible, if the generated number is less than 0 or bigger than 1, a new random number between 0 and 1 is generated and replaces with the infeasible one [44].

B. RKBMO FOR BAP
For BAP solution representation, the number of vessels is considered as the number of genes in the bird (solution), where every vessel represents a gene [43], [41]. In order to assign vessel to a berth, we generate, for each vessel, a random integer number in {1, 2, . . . , m} and add a uniform number from [0, 1]. In decoding, the integer part of any random key represents the berth assignment for that vessel, and the fraction parts represent a vessel sequence for each berth. Fig. 3 shows an example of how the BAP is represented using a random-key scheme. The example consists of a BAP instance with 10 vessels and three berths. In the encoding step, a random integer number is generated from [1], [3] for each vessel, and a uniform (0, 1) deviate is added to the number. Again, the integer part of each number represents the berth assignment for that vessel, which means that the same integer part represents the set of vessels at the same berth, and the fractional part represents the vessel sequence for the berth. To decode the solution, we follow the steps as presented in [55] and [56]: 1) The berth number of each vessel is obtained based on the integer part of the random key (see Fig. 3).
2) The random keys are sorted for each berth in ascending order to find the vessel sequences.
3) The index of each vessel is obtained to represent the final solution of BAP.
During the search process, some numbers might be outside the feasible solution space. We handle this limitation by doing the following: if the value of the gene is less than zero, we randomly generate a random number, rnd between (0, 1) and add it to the first berth number (i.e. 1+rnd). If the value of the gene is bigger than the number of berths, we randomly generate a random number in (0, 1) and subtract this number from the number of berths (i.e. 3rnd) [43].

V. THE PROPOSED DISCRETE BIRD MATING OPTIMIZER
In order to propose a discrete version of the BMO, it is necessary to understand the main components of BMO and how the algorithm works. BMO has three main operators which it can use to produce a new solution. These are [14]: self-recombination, two-parent recombination and multi-parent recombination. These operators are implemented using the five mating strategies: monogamous, promiscuous, polyandrous, polygynous and parthenogenetic, as discussed in Section 2. In addition, the BMO applies a mutation operator with probability after the following strategies [14]: monogamous, promiscuous, polyandrous and polygynous. Therefore, the mutation operator is one of the main operators in the BMO that needs to be mapped onto a discrete space. VOLUME 8, 2020 Before presenting the proposed discrete version of the BMO, the following subsections describe these three discrete operators in more detail.

A. SELF-RECOMBINATION IN DBMO
In the original BMO, the self-recombination operator generates a new solution by modifying the genes of a selected solution [14]. One of the most popular operators that can substitute the basic self-recombination is the basic version of the hill-climbing (HC) algorithm. The HC algorithm starts with an initial solution and improves it by iteratively generating a neighbor solution. The current solution is replaced by the neighbor solution if the quality of the neighbor solution is better than the current one. The HC algorithm stops when the termination criterion is met [13], [57]. Hence, the basic HC algorithm (see Algorithm 1) is similar to parthenogenetic improvement behavior, where the HC algorithm generates a new solution by making some changes to the genes of the current solution. However, parthenogenetic improvement requires a probability rate to decide either to perform the parthenogenetic improvement or to skip the procedure.

Algorithm 1 Pseudocode of the Modified HC Algorithm
Input: x as the starting solution, mcf p . Improve = true. Generate random number r between 0 and 1 If r < mcf p then Repeat Generate (N (x)). //generate a candidate neighbor// if x better than x then //x is the generated Output: final solution found x.

B. TWO-PARENT RECOMBINATION IN DBMO
In the original BMO, the two-parent recombination operator generates a new solution by mixing and recombining the genes of two parents [14]. This reproduction strategy is similar to the standard crossover operator where two parents produce a new child by combining their genes. In other population-based approaches, several crossover methods have been developed for solving combinatorial problems. For example, partially matched crossover [58], cycle crossover [58] and order crossover (OX) [58], [59]. Experimental results conducted by [60] revealed that OX performs better than others. Therefore, OX is used as a two-parent recombination operator in the DBMO. More details about the steps of OX are reported in [60].

C. MULTI-PARENT RECOMBINATION IN DBMO
The multi-parent recombination operator in the original BMO generates a new solution by combining the genes of more than two parents [14]. In order to convert the multi-parent recombination operator from a continuous space onto a discrete space, we need a multi-parent recombination operator that can deal with combinatorial problems. Two methods have been proposed in the literature to tackle this issue: multi-parent partially mapped crossover (MPPMX) [61] and adjacency-based crossover [62]. The experimental results presented in [62] demonstrated that adjacency-based crossover is not efficient in producing a good-quality solution for more than two-parent crossovers [62]. Therefore, the MPPMX was chosen in this work to substitute the multi-parent recombination operator in the original BMO. The MPPMX was proposed in order to extend the basic two-parent partially mapped crossover (PMX) into a multiparent crossover for better performance [61]. There are four main steps in MPPMX to construct a new offspring: substring selection, substring exchange, mapping list determination and offspring legalization [61]. The details of these steps are given in [61].

D. MUTATION OPERATOR IN DBMO
In the original BMO, the mutation operator is applied with probability after the following strategies: monogamous, promiscuous, polyandrous and polygynous [14] have been performed. The mutation operator is applied using an equation that can deal with continuous problems only (see Section 2). Therefore, an insertion operator is introduced in the DBMO to substitute the continuous mutation operator in the original BMO. In the insertion operator, an element is randomly selected to be removed and is added to another randomly selected position. Note that the same mutation probability (mcf) of the original BMO is applied in the DBMO.

E. THE DBMO FREMWORK
The following steps explain the DBMO procedure, which was adapted from the basic BMO procedures [14]: Step 1 (Parameter Initialization): DBMO uses the same parameters as the basic BMO (see Section II).
Step 2 (Society Initialization): randomly initialize a set of feasible solutions and add them to the society. The initialization strategy is related to the problem domain.
Step 3 (Society Evaluation): calculate the quality of each bird using the objective function of the problem.
Step 4 (Ranking): rank the birds in the society based on their quality.
Step 5 (Classification): classify the birds in the society into five groups based on their quality into: parthenogenetic, polyandrous, monogamous, polygynous, or promiscuous birds. The classification process is the same as in the original BMO (see step 5 in Section II).
Step 6 (Breeding): each bird produces a new brood using its own pattern, see Fig. 4. Step 7 (Replacement): If the quality of the brood is better than the quality of the bird, the brood replaces the bird. Otherwise, the bird remains in the society, and the brood is removed.
Step 8 (Termination Condition): repeat steps 4 to 7 repeated until a predetermined number of generations (NG) is performed.
Step 9 (Report the Best): select the bird with the best quality in the society as the best solution.
The pseudocode of the DBMO is illustrated in Algorithm 2.

VI. DBMO FOR TSP AND BAP
In this section, we present the solution representation for both BAP and TSP using the DBMO.

A. DBMO FOR THE TSP
In this study, the solution is represented using the simple popular representation for the TSP, i.e. a one-dimensional array with the length of N , where N is the number of cities as in [43,68]. Each gene in the solution represents a city number.  In this work, the TSP solution is generated in a random manner. First, we generate an array that contains the total number of cities. Next, we generate an empty array, randomly pick a number from the first array, and add it to the empty array. The selected number is then removed from the first array. This process is repeated until all the cities are added into the empty array.

B. DBMO FOR THE BAP
In the BAP, the solution is presented using a string of integer numbers as in [64]. The length of the solution is the number

Algorithm 2 Pseudocode of the DBMO Algorithm
Determine the society size (SS), maximum number of generation (), number of mates (nm) and mutation control factor (mcf). Generate SS feasible birds. for t = 1 to gen max do Rank the birds in ascending order based on their quality. Classify the society into five groups: parthenogenetic polyandrous, monogamous, polygynous and promiscuous. //the classification is based on the quality of each bird. of vessels n plus the number of berths (m) minus one. That is, the integer numbers contain m segments separated by 'zeros'. Each segment represents a service sequence for a number VOLUME 8, 2020 of particular vessels at the assigned berth. Fig. 6 shows an example of the solution representation of an instance of the BAP that contains 13 vessels and three berths. The solution has a length of 15 and contains two zeros to separate the berths. For example, vessels 4, 3, 1 and 8 will be serviced at berth 1 in that order. The solution is generated in three steps (as in [32]). First, we create an array with the number of vessels n and zeros m-1. Next, an empty solution array is created and then the elements of the first array are randomly selected and moved to the solution array. Finally, the vessel sequence for each berth is sorted in ascending order based on vessel arrival time. The pseudocode of the DBMO initialization steps for the BAP solution is given in Algorithm 3.

VII. COMPUTATIONAL RESULT
In this section, we discuss the results of an evaluation of the effectiveness of the proposed BMO variants RKBMO and DBMO and identify which of the two variants performs best. We also compare the performance of the proposed algorithms with that of the basic GA and PSO [65], [66]. The GA is chosen for this comparison because the proposed RKBMO and DBMO are evolutionary based and the GA is the standard algorithm for evolutionary algorithms [67], [68], while PSO has been widely improved and used to address many combinatorial optimization problems [24], [43], [46], [69], [70]. The GA has been implemented using roulette wheel selection with order crossover [71]. All algorithms (RKBMO, DBMO, GA and PSO) are tested on two combinatorial optimization problems: TSP and BAP. A statistical analysis is also conducted to verify the obtained results. First, however, the parameter settings and the experimental design for both the RKBMO and DBMO are presented.

A. PARAMETER SETTINGS AND EXPERIMENTAL DESIGN
The proposed algorithms were implemented using Java NetBeans IDE version 8.1 on a personal computer (Intel Pentium (R) Core i5 CPU at 3.40 GHz with 4 GB RAM), running a Windows 10 operating system (64-bit). The parameter settings of both the original BMO and DBMO for TSP and BAP are presented in Table 1. To ensure a fair comparison of the tested algorithms, we used the same parameter settings for both the basic BMO and DBMO. The values of SS and NG were obtained based on preliminary experiments as follows: 200 and 4000 for the TSP and 100 and 2000 for the BAP, see Tables A1 and A2 in Appendix A. The remaining parameters were fixed in line with [14] (see Table 1). Likewise, for the GA and PSO, for a fair comparison, the maximum number of generations and the population size were set as the same as those for the RKBMO and DBMO. All algorithms (GA, PSO, RKBMO and DBMO) were executed over 30 independent runs with different random seeds for all instances of the problem domains. The reason for executing the proposed algorithms on 30 runs was to obtain a good indication of algorithm consistency, and to enable a robust statistical analysis of algorithm performance [13]. In the subsequent tables, the results for each instance are presented as average (Avg.), average time (Avg Time), standard deviation (Std), and percentage deviation (Gap%) with respect to the quality of solution produced by the compared algorithms. The best results are highlighted in bold and the Gap is calculated as follows: where BCA show the best obtained from the compared algorithms and BKS is stand for best known solution in the literature.

B. RESULTS FOR THE TSP DATASET
The experimental results of the GA, PSO, RKBMO and DBMO on the TSP dataset are presented in Table 2. In addition, the computational time of the DBMO was much better than that of the GA, PSO and RKBMO in most instances (13,20 and 20 out of 20 instances, respectively), but it was sometimes more expensive than the GA especially for large-size instances. The reason for the high computational time of the DBMO in some instances is that it uses multi-parent partially mapped crossover and local search, which improves the solution quality at the expense of increased computation time. However, the computational time of the DBMO increased in only seven out of 20 instances and the overall performance of DBMO was superior to that of the GA, PSO and RKBMO. Similar to the results for the reported for the BAP, these results indicate the positive impact of the improvement operators of the DBMO on performance.
In order to verify whether the DBMO is significantly different from GA, PSO and RKBMO, a multiple comparison test was performed using Friedman and Quade statistical tests at a significance level of 0.05 [72]. (Note that a non-parametric test was performed based on normality tests and indicated that the data was not normally distributed.) The Friedman test is performed first, and if significant differences are found (the P-value of Friedman or Iman-Davenport statistic is less than the critical level 0.05), the Friedman and Quade tests are conducted to calculate the average rank of the compared algorithms to estimate the best performing one (the lower the better).
Then, we perform post-hoc procedures to obtain the adjusted p-values for each comparison between the control algorithm (the first ranking one) and the rest. (more details are reported in [72]). The p-value computed through Friedman and Quade tests (0.0000 and 0.0000) and the Iman-Davenport (0.000) are less than the critical level 0.05, therefore, we performed post-hoc procedures to detect the significant difference between all tested methods using Holm and Hochberg statistical tests.  Table 3 shows the average ranking (the lower the better) of Friedman and Quade statistical tests for the compared VOLUME 8, 2020 methods on TSP datasets. DBMO obtained the first rank followed by PSO, GA and RKBMO, in the given order. Therefore, DBMO will be the controlling method in Holm and Hochberg statistical tests. The Holm and Hochberg statistical test of Friedman and Quade tests [72] is conducted to detect the significant differences for each comparison between DBMO, GA, PSO and RKBMO. Table 4   The distribution of the results obtained from the GA, PSO, RKBMO and DBMO for four TSP instances (eil51, ch130, ch150 and kroA200) after running the program 30 times is presented in Figure A1 (see Appendix A). The graphs present the maximum and minimum values. The boxes represent the centre 50 percent of the data. The lower quartile of the boxes shows the lower median of the data while upper quartile presents the median of the upper part of the data. The line across the boxes relates to the median of the data. This is to show the distribution of the solutions obtained from the GA, PSO, RKBMO and DBMO. The figure shows that the distribution of the solutions is symmetric and distributed around the median in most instances. Based on these graphs, it can be seen that, for most instances, the DBMO is more consistent than the GA, PSO and RKBMO.

C. RESULTS FOR THE BAP DATASET
The results reported in Table 5 indicate that the DBMO, across all instances, outperforms the GA, PSO and RKBMO not only in terms of the best solution quality, but also with respect to the average quality and the average of all standard deviations of GA, PSO, RKBMO and DBMO is equal to 54.77, 36.20, 222.17 and 7.42. In addition, the DBMO finds the best-known solution (Gap% = 0) in seven instances (i03, i04, i05, i14, i15, i17 and i18) while the RKBMO and GA didn't find any best-known solution. (best-known solution here is the optimal value in the literature and can be found at [73]). This result demonstrates the effectiveness of the improvement process of the DBMO, which uses different discrete operators (two-parent, multi-parent and local search) to modify and improve the solution of the problem at hand. The RKBMO, on the other hand, uses mathematical equations to search the space and improve the solutions. The basic GA only uses two-parent crossover and mutation operators. The GA took less computational time than the DBMO and, basically because the use of multi-parent crossover in DBMO algorithm consumes extra computational time.
We performed a multiple comparison test between GA, PSO, RKBMO, and DBMO using Friedman and Quade statistical tests. Table 6 shows the average ranking (the lower the better) computed by Friedman and Quade tests for all compared methods. The table highlighted that DBMO obtained the first rank out of four compared methods followed by PSO, GA and RKBMO, respectively. The p-value computed through Friedman and Quade tests are less than 0.05 (0.0000 and 0.0000) and the Iman-Davenport (<0.05) prove that there is a significant difference among the methods. Therefore, post-hoc procedures (Holm and Hochberg) are performed to detect the significant difference between all tested methods. The adjusted p-values of Holm and Hochberg statistical tests presented in Table 7 (where DBMO is the controlling method) demonstrate that DBMO outperforms all other methods with critical level of 0.05 (adjusted p-values < 0.05).
The box-whisker plot in Fig. A2 (see Appendix A) shows the distribution of the results for 30 runs of the DBMO, RKBMO and GA for four instances of the BAP (i01, i02, i03 and i04). It can be seen from these graphs that the distribution of the results for most of the instances in the case of the DBMO is skewed to the lower end (i.e. i01, i03 i04). However, there are some cases where the solution distribution is symmetric; the solutions are distributed evenly around the median (i.e., i02). While the distribution of the solutions obtained by the GA, PSO and RKBMO in most instances is symmetric, the distribution is split at the median. In addition, in the case of the RKBMO one solution distribution is skewed to the upper end (worst solution) (i01). These distributions demonstrate that the DBMO is more consistent than the GA, PSO and RKBMO.
To conclude, it is clear that in both problem domains (TSP and BAP), the DBMO outperformed the GA, PSO and RKBMO and matched the best-known results in some instances. In addition, for both domains the standard deviation was relatively small. Also, the percentage deviation (Gap) showed that the DBMO results were very close to the best-known solutions in both domains. This favorable outcome reveals the benefit of using the discrete version of the BMO rather than the original one with random-key representation. The DBMO was able to achieve a significant difference in performance compared to the RKBMO for all instances of both the TSP and BAP. This is due to the improvement procedures in the DBMO: two parent recombination using OX and multi-parent recombination using MPPMX which gives a possibility to combine the components of many solutions and generate new one. Another feature of DBMO is the use of Hill Climbing heuristic for self-recombination which improves the ability 96854 VOLUME 8, 2020   of the DBMO for exploiting the most promising regions of the search space. In contrast to the RKBMO, which uses mathematical equations through random-key representation of the solution to move in the search space and improve the solution, the DBMO uses discrete operators that can directly deal with the solution and improve it. The superior results produced by the DBMO clearly answer the question of which variant of the BMO that should be used for combinatorial optimization. problems as well as the best operators and configuration for applying it to such problems.

VIII. CONCLUSION
In this work, we proposed two variants of the BMO for solving combinatorial optimization problems: the random-key bird mating optimizer (RKBMO) and discrete bird mating optimizer (DBMO). The RKBMO used random-key representation to construct a direct mapping relationship between the algorithm and the problem at hand. In contrast, the DBMO used three main operators to generate new solution: the hill-climbing algorithm, two-parent order crossover, and multi-parent order crossover. To evaluate the proposed methods, we tested them on the travelling salesman problem and berth allocation problem using the same parameter settings. In addition, a comparative study was conducted with two widely studied optimization algorithms in the literature: GA and PSO. The computational results demonstrated that the DBMO was able to outperform the RKBMO across all instances in both problem domains as well as GA and PSO. Furthermore, the DBMO matched the best-known results in some instances. In conclusion, we proved that the DBMO can be used to solve combinatorial optimization problems and achieve satisfactory results. VOLUME 8, 2020 In future work, it might be beneficial to investigate the effectiveness of combining the three main DBMO operators (two-parent crossover, multi-parent crossover, and local search) in a pool in order to select one of them randomly and employ it to improve the whole population. In addition, changing the selected operator to another one could be done according to the search status. This may improve the DBMO to allow it to effectively tackle any problem with minimal modification.
ANAS ARRAM received the Ph.D. degree in computer science with The National University of Malaysia, in 2017. He is currently working as a Computer Vision Engineer with FitLens Company. He has published more than six papers at international journals peer-reviewed international conferences. His main research areas include meta-heuristics, hyper-heuristics, scheduling and timetabling, machine learning, and deep learning. He is a member of the Data Mining and Optimization Research Group (DMO), Centre for Artificial Intelligent (CAIT), National University of Malaysia. He received an Outstanding Researcher Award upon completing his Ph.D.
MASRI AYOB received the Ph.D. degree in computer science from the University of Nottingham, in 2005. She was a member of the ASAP Research Group, University of Nottingham. She has been a Lecturer with the Faculty of Information Science and Technology, Universiti Kebangsaan Malaysia (UKM), since 1997. She is currently a Principle Researcher with the Data Mining and Optimisation Research Group (DMO), Centre for Artificial Intelligent (CAIT), UKM. She has published more than 150 articles at international journals and at peer-reviewed international conferences. Her main research areas include meta-heuristics, hyperheuristics, scheduling and timetabling, especially educational timetabling, healthcare scheduling, and routing problems. She has been served as a program committee for more than 50 international conferences and a Reviewer for high impact journals.
GRAHAM KENDALL received the B.Sc. degree (Hons.) in computation from the University of Manchester Institute of Science and Technology (UMIST), U.K., in 1997, and the Ph.D. degree from the School of Computer Science, University of Nottingham, in 2000. Before entering academia, he spent almost 20 years in the IT industry, working for various U.K. companies (including the Co-operative Wholesale Society and Provincial Insurance), undertaking a variety of roles including a Computer Operator, a Technical Support Manager, and a Service Manager. He is currently a Professor of computer science. He is also with the University of Nottingham's Malaysia campus, where he is also serving as the Vive-Provost of research and knowledge transfer. He is also a member of the Automated Scheduling, Optimisation and Planning Research Group, School of Computer Science. He has edited and authored 9 books and has published almost 50 refereed journal articles (the vast majority in ISI ranked journals) and more than 90 articles in peer reviewed conferences. His expertise includes operational research, meta-and hyper-heuristics, evolutionary computation, and artificial intelligence, with a specific interest in scheduling, including timetabling, sports scheduling, cutting and packing, and rostering. He is a Fellow of the Operational Research Society. He has been awarded externally funded grants worth over £5.5M from a variety of sources including the Engineering and Physical Sciences Research Council and commercial organizations. He chairs the Steering Committee of the Multidisciplinary International Conference on Scheduling: Theory and Applications, in addition to chairing several other international conferences in recent years. He is an Associate Editor of eight international journals (including two the IEEE TRANSACTIONS).
AlAA SULAIMAN received the B.A. degree in electrical and computer engineering from AN-Najah National University, Nablus, Palestine, in 2003, and the M.Sc. degree in scientific computing minor in computer science and engineering from Bierzit University, Ramallah, Palestine, in 2010. He is currently pursuing the Ph.D. degree with the Faculty of Information Science and Technology, Universiti Kebangsaan Malaysia (UKM). His area of interests include machine learning, image processing, and pattern recognition.