Hybrid Bird Mating Optimizer With Single-Based Algorithms for Combinatorial Optimization Problems

Bird mating optimizer (BMO) is a population-based metaheuristic that has been recently extended to solve combinatorial optimization problems. Even though the algorithm shows promising performance in solving combinatorial optimization problems, it suffers from slow convergence and poor efficiency which leads to poor solution quality for some problem instances. This is due to the limited capability of BMO in exploiting the search space and identifying more promising regions. Therefore, in this work we propose a hybrid BMO with five single-based metaheuristics: hill-climbing, late acceptance hill-climbing, simulated annealing, iterated greedy heuristic and variable iterated greedy heuristic. Each of these algorithms is used inside the BMO to exploit the search space, and improve the quality of solution generated from the BMO population. This work also compares which one of these five is better for hybridizing with BMO. The performance of these algorithms is tested on two combinatorial problems: travelling salesman problem and berth allocation problem. Experimental results demonstrate that the hybrid algorithm is superior to BMO when applied to both problems and it improved the BMO by 1.13% for BAP and by 4.13% for TSP. Furthermore, the hybrid algorithm is able to match the best-known results for most of the instances. In addition, the proposed hybrid approaches perform well over both tested domains and obtain competitive results when compared to the best-known results that have previously been presented in the scientific literature.


I. INTRODUCTION
Combinatorial optimization problems (COPs) arise in many areas such as computer since, operational research, electronic commerce, and artificial intelligence. COPs can be defined as a topic to find an optimal ordering or arrangement for given discrete variables [1], [2]. Examples of these problems are job scheduling, journey planning, university educational timetabling, travelling salesman, and berth allocation [3]. COPs are NP-hard which may not be suited for exact methods that aim to find the actual optimal solution. Instead, approximate algorithms that can find (near) optimal solution with reasonable time are more practical especially when the problem instance is sizable [4]. Approximate algorithms consist The associate editor coordinating the review of this manuscript and approving it for publication was Jamshid Aghaei . of heuristic and metaheuristic algorithms in which heuristic is a problem-dependent algorithm that is designed to solve a specific problem, while metaheuristic can work in effectively across different problems. Therefore, metaheuristics are the better choice when solving several COPs. Examples of these algorithms include simulated annealing (SA) [5], hillclimbing (HC) [6], late acceptance hill-climbing (LAHC) [7], genetic algorithm (GA) [8], [23], harmony search [9], iterated greedy algorithm (IG) [10]- [12], ant colony optimization [13], [14] and particle swarm optimization [15].
Bird Mating Optimizer (BMO) is a natural inspired metaheuristic algorithm that was proposed by [16]. The BMO inspired by mating behavior of birds during mating season, in which each male bird attempts to mate with a female one to breed a new brood. BMO has some advantages as avoiding trap in local optima by using five mating strategies to move through the search space. In addition, BMO has more capabilities in exploring and exploiting the search space than other mating algorithms such as genetic algorithm [16]. BMO shows a competitive result comparing to other evolutionary algorithms such as classical evolutionary programming, fast evolutionary programming, classical evolutionary strategies, fast evolutionary strategies, genetic algorithm, particle swarm optimization, and group search optimizer [16]. However, the BMO was originally proposed to solve continuous optimization problems, but in 2020 [17] proposed a discrete version of BMO for solving COPs. In this work, we will use this version of BMO as it shows an effective and promising performance for solving COPs.
However, even though BMO shows promising performance in solving optimization problems, it suffers from slow convergence and poor efficiency for some constraint problems. This is due to the limited capability of BMO in exploiting the search space and identifying more promising regions. For some complicated problems, BMO shows premature convergence or poor efficiency [18], [19]. Therefore, maintain the balance between exploration and exploitation ability of the BMO algorithm is needed. Some researchers have attempted to tackle these issues by conducting several improvement modifications. One of these modifications is hybridizing BMO with single-based metaheuristics [19], [18] to improve the exploitation ability of the algorithm. In this study, our aim is to hybrid the BMO with five singlebased metaheuristics, HC, LAHC, SA, IG and variable IG, to improve the performance of BMO for solving COPs. The reason for choosing five single-based metaheuristics is to compare and investigate which of these five can lead to better performance when hybridizing it with BMO.
According to Talbi et al. [6], metaheuristic hybridization can be done by combining a metaheuristic with a complementary metaheuristic. In this model, there are two levels of hybridization: low-level and high-level. In the low-level hybridization, a given function (e.g., crossover or mutation in GA) of metaheuristic is replaced by another metaheuristic (e.g., single-based metaheuristic). While in the high-level, the different metaheuristics are self-contained and there is no relationship between their internal working. For each level, two hybridization mechanisms are available: relay and teamwork hybridization. In relay hybridization, a set of metaheuristics is applied one by one and each one is taking the output of the previous one as input. Whilst, teamwork hybridization represents a set of optimization models that cooperated as agents in parallel, and each agent searches in solution space [6].
In this work, we aim to use the low-level teamwork hybrid (LTH) model for COPs, where a population-based algorithm BMO will be hybridized with five single-based algorithms, respectively. In other words, each single-based algorithm will be improving the exploitation ability in BMO and identifying more promising solution regions. The single-based algorithm will be applied with a probability to every newly generated solution in the population of BMO to ensure that it's the local optima. In addition, a comparative study will be conducted between these five hybrid models to figure out which one is more effective in finding more promising solutions.
The rest of the paper is organized as follows: Section 2 represents the related work. The basics of BMO, HC, LAHC, SA, IG, and VIG are presented in Section 3. Section 4 presents the details of the proposed hybridization approach. Section 5 summarizes and analyzes the experimental. Finally, the conclusion and details of future works are given in Section 6.

II. METHODS
The following subsections present the details of the BMO algorithm and the five single-based algorithms.

A. BIRD MATING OPTIMIZER
BMO is a population-based metaheuristic algorithm proposed by Askarzadeh [16] that mimics the mating behavior of bird species during mating season. The population is referred as a society and each member of the society represents a feasible solution and called bird. The new generated solution called brood. The society contains of two components: males and females. Males are classified into three components: monogamous, polygynous, and polyandrous. While females are with the most promising genes which contain two categories: parthenogenesis and polyandrous. In total, the algorithm uses five updating strategies to generate new solutions and they are explained below with details. Note here that we are using the BMO version that proposed by Arram et al. [20] which was proposed to solve COPs.
Parthenogenesis is the mating system in which the female bird can produce a brood without mating with a male. In this system, 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 using simple non-improvement hill-climbing (NI-HC). The NI-HC algorithm starts with a solution and generates a new one by applying some changes to the genes and repeat this process until there is no improvement to the current solution.
Applying this system requires a probability at every iteration to decide to apply this algorithm or to skip it. The pseudocode of this system is presented in Algorithm 1 [16], [17].  Monogamy is a mating system in which a male tends to mate with only one female. Every male selects its mates from the female group (parthenogenesis and polyandrous) by evaluating the quality of the females using a probabilistic approach to select one of them. So female birds with good genes have a higher probability of being selected. For this system, the two-parent order crossover (OX) is used as a mating system between the male bird and its selected female [16], [17]. More details about OX are reported in [21]. The produced brood (new solution) from this mating system, will go through a mutation operator with a probability to improve the quality of the brood. Insertion operator will be used as a mutation operator where a gene is selected from the brood randomly and moved to another randomly selected position in the brood.
Polygamy is a mating system in which each polygynous male 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 birds mate with several females to produce 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 multi-parent partially mapped crossover (MPPMX) is used in this system as a mating operator between males and several selected females. MPPMX was proposed by Ting et al. [22] as an extension of two-parent partially mapped crossover to multiparent crossover for better performance. The MPPMX applies four main steps to generate a new brood: substring selection, substring exchange, mapping list determination and offspring legalization. More details of this crossover are reported in [22]. In addition, insertion mutation operator is applied to the produced brood with a probability to improve its quality [17].
Promiscuity is also another mating system in which one male mates with several females with an unstable relationships. This mating process indicates a chaotic social structure in which the male bird will never see the brood, 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. Therefore, the OX is used in this mating process and insertion mutation operator will be applied with a probability [17].
Polyandry is the last mating system where a polyandrous female bird seeks to mate with more than two monogamous males. The female performs a selection mechanism to select the males. Then, each female bird mates with her selected male birds. MPPMX is used for this mating system as multiple birds will mate, and insertion mutation operator will be applied with a probability [17]. Note here that the mutation operator with a probability is applied to the four mating systems (Monogamy, Polygamy, Promiscuity and Polyandry) regardless of the resulting solution quality, and the aim of this is to maintain the diversity of the algorithm.
The following steps explain the BMO algorithm procedure as presented in [17]: Step 1: Parameter initialization: initialize the BMO parameters. Step 2: Society initialization: initialize the population randomly with feasible birds (solutions). 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 in descending order.
Step 5: Classification: classify the birds in the society into five groups based on their quality, the classification is as follow: parthenogenetic, polyandrous, monogamous, VOLUME 9, 2021 polygynous, and promiscuous birds. The percentage of each group is a parameter and defined in the parameters initialization step.
Step 6: Breeding: each bird produces a new brood using its own mating system.
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 society, and the brood is removed.
Step 8: Termination condition: steps from 4 to 7 are repeated until a predetermined number of generations is performed.
Step 9: Report the best: select the bird with the best quality in society as the best solution.
The pseudocode of the DBMO is illustrated in Algorithm 2.
Insert the current cost into the fitness array Hill Climbing (HC) is one of the well-known local search methods that attempt to find better solution. The algorithm starts with an initial solution and improves it by generating a neighborhood solution. The current solution is replaced with a neighborhood solution if the quality of the neighborhood is better than the current, otherwise, the neighborhood is rejected, and HC begins a new iteration. The search process continues until the stopping criteria is satisfied [23], [24]. Algorithm 3 presents the pseudo-code of the algorithm.

C. LATE ACCEPTANCE HILL-CLIMBING
The late acceptance hill-climbing algorithm (LAHC) is an improved version of HC algorithm, which was proposed by Burke and Bykov [25]. The main idea of LAHC is to accept the non-improving solution when the quality of the new generated solution is better (or equal) than those that were recently accepted a few iterations before. In practice, LAHC starts from a given initial solution, and iteratively improves it by comparing the new candidate solution with the current one in order to accept or reject it. Hence, to apply the LAHC rule, the algorithm will create a list with a fixed length to

Algorithm 3: Pseudo-Code of SA Algorithm
Generate the initial solution x = x 0 Set the starting temperature T = T 0 Set Cooling rate β. Set max number of trials at each temperature i max Generate a random number between 0 and 1 (r) save the quality of recently visited solutions. If the quality of the new candidate solution is better than the quality value of the last element in the list, the candidate solution will be accepted as the current initial solution. Then, the element at the end of the list will be removed and the quality value of the newly accepted solution will be added to the beginning of the list [7]. Algorithm 4 presents the pseudo-code of the LAHC algorithm [7], [25].

D. SIMULATED ANNEALING
Simulates annealing (SA) was proposed by Kirkpatrick et al. [5], Černý [26] as an optimization algorithm for solving optimization problems. SA was proposed based on hill-climbing to overcome its problem of trapping in local optima. SA uses probability to accept worse solution and give more chances to explore the search space. The algorithm starts with a random initial solution, and at each iteration, a neighbor solution is generated using a predefined neighborhood structure. The new solution is evaluated using fitness function and will be accepted if it is better than the best one, otherwise, the worse solution is accepted with a probability that is determined by Boltzmann probability as in equation 1. In addition, T is a temperature that periodically decreases during the search process according to the cooling schedule [5], [27]. The pseudocode of SA is presented in Algorithm 5.
where E is the difference between the best solution and the generated one.

Algorithm 4: Pseudo-Code of IG algorithm
Generate the initial solution x Set the starting temperature T = T 0 Repeat //destruction procedure for i = 1 to d do //d is the destruction size remove one gene at random from x and add it in x2 // x2 is a second array to store the removed genes from x end //construction procedure for i = 1 to d do x = best permutation obtained by inserting gene x2 i in all possible position of x End if the quality of x is better than the x do Update the best solution Until Stopping criteria satisfied Return best found solution.

E. ITERATED GREEDY ALGORITHM
Iterated greedy algorithm (IG) is a stochastic search method that was proposed by [10]. The IG generates a new solution using the idea of destruction and construction phases. In destruction phase, number of solution's genes d are randomly chosen and removed from the solution without repetition, so that two partial solutions will be resulted. The first, with the size d of genes, is donated as SR including the removed genes in the same order that they removed. The second, with the size n-d of genes, is the original solution without the removed genes which is donated as SD. Next, the construction phase is employed to reinsert the removed genes into the solution. The NEH insertion heuristic is used as a constructive procedure to complete the solution. The first gene SR1 is inserted into all possible n -d+ 1 positions in the destructed solution SD, which generates n -d+ 1 partial solutions. Among these n -d+ 1 generated partial solution, only the solution with best quality is chosen and kept for the next iteration. The second gene is then considered, and the process continues until SR is empty or final solution is obtained. Therefore, SD is again with the size of n [10]. In addition, the Boltzmann probability of SA is used here (see equation 1), where new generated solution is always accepted if it is better than the best, otherwise the worse one is accepted with the probability. But one different thing that is the temperature Update the best solution Until Stopping criteria satisfied Return best found solution.
here is static, means there is no cooling schedule [10]. The pseudocode of IG is given in Algorithm 6 [10].

F. VARIABLE ITERATED GREEDY ALGORITHM
Variable iterated greedy algorithm (VIG) was proposed by Framinan and Leisten [28] to solve the permutation flowshop scheduling problem. The VIG is inspired from the variable neighborhood search (VNS) algorithm which presented in [29]. The VIG is developed by using the idea of neighborhood change of the VNS algorithm. The implementation of VIG is similar to IG but differs in the destruction parameter, where the maximum destruction size is fixed at d max = n-1. The destruction size is initially set to d = 1. The current solution is destructed and reconstructed again with the variable size of d. Then, destruction size d is incremented by 1 (i.e., d = d+1), if the solution is not improved until d max = n-1. Whenever the solution is improved in any destruction size, the destruction size is again set to d = 1 and the search starts from the beginning once again. The pseud-code of the VIG algorithm is given below [28]. VOLUME 9, 2021 Algorithm 6: 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 < probability then Repeat Generate (N (x)). //generate a candidate neighbor// if x better than x then //x is the generated Output: final solution found x.

III. PROBLEMS DESCRIPTION
In this work, the performance of the proposed hybrid BMO 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 is classified as an NP-hard problem, which means that it may take an infeasible computational time to solve it [30]- [32]. The TSP can be described as a search for the shortest path route between a list of cities, that visits each city once and only once, and finally return to the origin city [33].

B. BERTH ALLOCATION PROBLEM
BAP is a NP-hard problem and known as a berth scheduling problem, where it aims to allocate berth space for vessels in container terminal to be served [34], [35]. The BAP can be classified into two models: the first one is discrete and continuous berth space, whilst the second is dynamic and static vessel arrival time. The berth space is considered as discrete if the quay is divided into a set of segments (berths), and continuous if not partitioned and vessels can berth anywhere along the quay. The vessel arrival time is considered as dynamic if vessels can arrive at any time during the container operations with planning arrival time, and static if all the vessels have to arrive in the port before the berth planning step begins. In this study, we focus on the discrete and dynamic version of the BAP [36]. This BAP deals with allocating vessels to berths in the port 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 [37]: 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 [44]: where: a i : arrival time of vessel i -K : set of berths, -N : number of vessels that will arrive in the harbour, -P k i : handling time of vessel i at berth k -T k i : berthing time of vessel i at berth k x k ij : decision variable, x k ij = 1 if vessel j is serviced at berth k immediately after vessel i

IV. THE PROPOSED APPROACH
In the hybridization process, invoking single-based algorithms to the whole society (population) of BMO would be computationally expensive and also might cause a premature convergence in the BMO in the early stage of the search [38]. Therefore, the probability of hybridizing every member in the society is controlled based on a probability rate (hybrid probability (HP)) which is tuned based on preliminary experiments. As mentioned in Section 1, we are using a low-level teamwork hybrid, therefore the hybridization process is as follows. At each iteration, BMO generates a new society using its five mating systems (monogamous, polygynous, polyandrous, parthenogenesis, and polyandrous). Next, all birds in this new society will go through a hybrid probability rate, and any bird passes the probability rate, will go through a single-based algorithm, otherwise the bird will remain the same in the society. Finally, the society will go to society update step in BMO. The flowchart of the proposed hybrid algorithm is presented in Figure 1. The Figure 1 shows the process of hybridizing BMO with HC, however, the same process is applied to the remaining four single-based algorithms (LAHC, SA, IG, VIG).
Furthermore, the new solution that obtained from the single-based heuristic algorithms is generated as follows. For IG and VIG, the new solution is generated based on the greedy manner as explained in sections 3.4 and 3.5, this is because they are designed to generate a new solution based on that greedy manner. Whilst, the new solution generated from HC, LAHC and SA is generated using neighborhood structure. Hence, the type of the neighborhood structure depends on the problem domain. For BAP, the new solution is generated using swap operator, where two vessels at any berth are randomly selected, and their positions are exchanged [36]. For TSP, the neighborhood structure is generated using 2-opt neighborhood operator. The 2-opt operator selects two tours at random from the current solution and swaps the cities located at the end sections of the selected tours. This operator is different from the normal swap in the way that the swap operator only exchanges the places of the two selected cities.

V. EXPERIMENTS
The following sub-sections present the experimental setups, parameters tuning, and results.

VI. DATASETS
In order to assess the performance of the proposed BMO variants, the experiments are performed over two COPs: TSP and BAP. TSP benchmark includes 20 instances that range from 51 to 318 cities [13]. Whilst BAP involves 30 different instances, each instance has 30 vessels and 13 berths [35].

VII. PARAMETER SETTINGS
The proposed algorithms were implemented using Java Net-Beans 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). All algorithms were executed over 30 independent runs with different random seeds for all instances of the problem domains.
The parameter settings of the BMO and single-based algorithms are set s follows:

A. BMO PARAMETER SETTINGS
Due to the impact of hybridizing single-based algorithms with BMO on its performance and diversity, the society size is carefully tuned. Therefore, the population size (PS) of the hybrid DBMO is fixed to 30 based on preliminary experiments for both TSP and BAP. The number of generations (NG) is also tuned based on preliminary experiments and set to 4000 and 2000 for TSP and BAP, respectively (see Table 10, Appendix A). The remaining parameters of BMO are fixed as suggested in [16] and presented in Table 1.

B. SINGLE-BASED PARAMETER SETTINGS
In this experiment, the stopping conditions of all single-based algorithms for both BAP and TSP are fixed based on the exe-cution time. Based on the preliminary test, the execution time is fixed to 0.2 sec for BAP and to 0.0015 sec * number of cities for TSP. Therefore, the results of the proposed algorithms can be fairly compared with each other. For the remaining parameters of the single-based algorithms, the settings are as follows: • For HC parameters, HC has only one parameter, the maximum number of iterations. This parameter has been set based on the execution time as mentioned.
• LAHC has a list size (L) parameter, which is fixed to 10 for BAP and 50 for TSP based on preliminary experiments.
• SA has another two parameters in this work: Initial temperature (T ) and cooling rate (β). The β is fixed at 0.85 and the T is fixed to be 50% of the value of the initial solution as suggested in [39], [40] for both BAP and TSP. Note that that loop will stop when the T reaches 0.
• IG has also a destruction size parameter (d), which is randomly selected in range of 4 to 7 as suggested in [41] for both TSP and BAP domains.
• VIG has another parameter which is the temperature level T . This parameter is fixed to 1000 for TSP based on preliminary tests, and to 0.05 × TP/NBS for BAP as suggested in [41]. Where TP denotes total handling time (service time of a vessel on a berth), which is the summation of handling times of all vessels; NBS is the total number of allowable vessels for each berth, and is the number of vessels. For each berth, if the vessel can be served, it is considered an allowable vessel for the berth.

VIII. RESULTS AND DISCUSSION
All results obtained from the proposed hybrid BMO variants are presented in this section. The basic BMO are compared to it hybridized variants to assess the effect of hybridizing single-based algorithms to the basic BMO. To find out the best hybridized approach among the proposed five approaches, all approaches are compared together in one table. In addition, the best approach is compared to the state-of-the-art algorithms to verify the effectiveness of the proposed approach. Statistical analysis is also conducted to verify the significance of the best hybridized BMO variant over other variants. In the subsequent tables, the results for each instance are presented as average (Avg.), standard deviation (Std), the best (Best), 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 stands for best-known solution in the literature. In this section, we investigate the effectiveness of the hybrid BMO with single-based algorithms by comparing the results of the basic BMO with its hybrid approaches. The hybrid approaches with BMO are as follows: BMO-HC, BMO, LAHC, BMO-SA, BMO-IG, and BMO-VIG. In addition, we compare the five hybrid BMO approaches against each other to investigate which hybridized algorithm is better for each problem domain. Table 2 and Table 11 (see Appendix A) present the results of BMO vs. BMO-HC, BMO-LAHC, BMO-SA, BMO-IG and BMO-VIG for the BAP. Whilst, Table 3 and Table 12 (see Appendix A) report the comparison results of BMO vs. BMO-HC, BMO-LAHC, BMO-SA, MO-IG and BMO-VIG for the TSP. Best results are shown in bold.
From Tables 1, 2 and 11 (see Appendix A), we can make the following observations: For BAP problem, across all instances, all hybrid BMO algorithms outperformed BMO not only in terms of percentage deviation, but also on average, standard deviation (see Table 2 and Table 10). We can observe that hybrid BMO with single-based algorithms enhanced the exploitation of BMO and produced better results. In comparison between all proposed hybrid BMO approaches, it's clear that BMO-HC performed better than BMO-LAHC, BMO-SA, BMO-IG and BMO-VIG not only in terms of percentage deviation, but also, on average results. Table 11 (see Appendix A) also shows that the standard deviation produced by BMO-HC is better than BMO-LAHC, BMO-SA, BMO-IG and BMO-VIG on 28 out of 30 instances. In addition, BMO-HC improved the BMO by 1.13% on the average of all instances in terms of solution quality (as shown in the last column in Table 2). On average computational time as can be seen in Table 12 (Appendix A), the values obtained from the BMO are less than those obtained from the hybrid BMOs. This is because the BMO suffers from poor exploitation ability, which led the algorithm to converge at an early stage of the search space. Thus, the best results are obtained with less time but with poor quality solution. In addition, the hybrid DBMOs require more time due to the use of SBH algorithms during the search, and this is still acceptable as the hybrid algorithms still run in a reasonable amount of time. Meta-heuristic techniques aim to provide a good quality solution with a reasonable amount of time, and this is what the hybrid BMO in this paper has achieved. For TSP, the results reported in Table 3 and Table 12 demonstrated that BMO-HC, BMO-LAHC and BMO-SA obtained better results for all tested instances, in terms of the percentage deviation, average, standard deviation. Whilst, BMO-IG and BMO-VIG outperformed BMO on 17 and 13 out of 20 instances in terms of average results, respectively, and on 9 and 8 out of 20 instances in terms of standard deviation, respectively. Comparing the hybrid approaches against each other, we can see that BMO-SA outperformed BMO-HC, BMO-LAHC, BMO-IG and BMO-VIG on 19 out of 20 instances. On average and results, BMO-SA outperformed the four compared hybrid BMO approaches in 10 out of 20 instances (see Table 3). Table 13 (see Appendix A) also shows that BMO-SA produced better standard deviation compared to BMO-HC, BMO-LAHC, BMO-IG and BMO-VIG on only 5 out of 20 instances. In addition, BMO-SA improved the BMO by 4.13% on the average of all instances in terms of solution quality (as shown in the last column in Table 3). For computational time, as shown in Table 13, all hybrid BMO algorithms consumed more computational time compared to BMO. That is due to the early convergence in BMO and to the use of SBH in the hybrid approaches.
To verify whether the results are statistically different, we have conducted a pairwise comparison using Wilcoxon test with significant level of 0.05 for comparing BMO against hybrid approaches (Tables 15 and 15)7 and a multiple comparison statistical test to compare the hybrid approaches to each other (Tables 4-7). Table 15 and Table 16 show the pvalues computed for the tested instances of the BAP and TSP, respectively. In this table, ''+'' indicates hybrid DBMO is statistically better than BMO (p-value < 0.05), ''-'' indicates hybrid DBMO outperformed by DBMO (p-value > 0.05) and ''='' indicates both hybrid BMO and BMO have the same performance (p-value = 0.05).
Based on p-value reported in Table 15 and Table 16, it can be concluded that BMO-HC, BMO-LAHC, BMO-SA, MO-IG and BMO-VIG results are statistically significant and better than from those produced by BMO across all instances of BAP problem domain (see Table 15). For the TSP, the pvalues demonstrate that BMO-HC, BMO-LAHC and BMO-SA are statistically better than BMO on all tested instances, whilst BMO-IG and BMO-VIG are statistically better than BMO on 16 and 12 out of 6 instances, respectively (see Table 16). VOLUME 9, 2021    For multiple comparison between the hybrid approaches, Table 4 and Table 5 show the average ranking (the lower the better) computed by Friedman, Friedman Aligned and Quade tests for the hybrid BMO approaches (BMO-HC, BMO-LAHC, BMO-SA, MO-IG and BMO-VIG) for both BAP and TSP domains, respectively. The tables highlighted that BMO-HC and BMO-SA obtained the first rank out of five compared methods for BAP and TSP problem domains, respectively. The p-value computed through the statistics of each of the test considered are less than 0.05 and the Iman-Davenport  (<0.05) prove that there is a significant difference among the methods considered for both BAP and TSP. Therefore, post-hoc procedures are performed to detect the significant difference between all tested methods. Table 6 shows the adjusted p-values of Holm and Hochberg statistical tests using the ranks computed by Friedman, Friedman Aligned and Quade tests, respectively where DBMO-HC is the controlling method.  Comparison between the proposed approaches and the state-of-the-art approaches for BAP in terms of best obtained solution (minimum cost).

TABLE 9.
Comparison between the proposed approaches and the state-of-the-art approaches for TSP in terms of best obtained solution (minimum cost).
considered for BAP. Whilst for TSP, it can be seen in Table 7 that the BMO-SA outperforms BMO-IG and BMO-VIG with a critical level of 0.05 (adjusted p-values < 0.05). However, the results in Table 7 indicate that BMO-SA does not outperform BMO-SA and BMO-LAHC (adjusted p-values > 0.10).
From the above comparison, we can conclude that hybridizing single-based algorithms with BMO enhances the VOLUME 9, 2021  search process and the performance of BMO. It can be stated that BMO with single-based algorithm is a robust algorithm that balances between exploration and exploitation when searching for global optimum. The single-based algorithm improves the exploitation ability of BMO and at the same time, it does not damage its exploration ability.
To show the hybrid algorithm improves the exploitation ability of the algorithm, we have conducted an analysis using one BAP instance as shown in Figure 2. The figure shows that  the hybrid BMO-HC maintained the exploitation ability of the algorithm and continuously improve BMO until the end of the search. Parts (a), (b), (c) and (d) in Figure 2 show the ability of the algorithm in exploring the search space followed by the ability of exploiting it (Figure 2 (e), (f), (g) and (h)). This behavior directs the algorithm toward a superior convergence.   In this figure, it can be clearly seen that the algorithm tries to explore the search space until 30% of the search space and then starts to exploit it after the 30% of the search space and continues the exploitation until finding the best solution and final convergence.
To summarize, the results achieved by hybrid BMO approaches show that these approaches perform better than the BMO. In addition, the hybrid approaches matched the best-known results in the literature on 30 out of 30 instances for BAP and on 17 out of 20 instances for TSP. From these results, it is clear that embedding single-based metaheuristics with the BMO outperform the BMO algorithm alone. In fact, that is mainly due to the ability of single-based meta-heuristics to improve the exploitation process of the algorithm and produce good quality solutions. In addition, the statistical tests revealed that BMO-HC outperforms the other hybrid BMO approaches for BAP, whilst BMO-SA outperforms BMO-IG and BMO-VIG for TSP. However, although the results of Holm and Hochberg suggest that BMO-SA is not statistically better than BMO-HC and BMO-LAHC for TSP nevertheless, the best and the average results shown in Table 3 and Table 11 reveals that BMO-SA produced better results compared to BMO-HC and BMO-LAHC.

B. COMPARISON WITH THE STATE-OF-THE-ART APPROACHES
In the previous section, we found from analyzing the results that BMO-HC and BMO-SA outperformed other hybrid approaches in terms of solution quality for BAP and TSP, respectively, on all tested instances. In this section, we compare the performance of the of the best hybrid approaches against the most related approaches in the BAP and TSP literature. The following subsections present the comparison for BAP and TSP. Table 8 presents the comparison results of percentage deviation according to solution quality between BMO-HC and the following approaches: Adaptive Large Neighborhood Search (ALNS) [42], Column Generation (CG) [43], Simulated Annealing with Restart Strategy (SA r.s ) [44], Tabu Search Algorithm (TS) [35], Iterated Greedy Heuristic (IG) [41]. Best results are highlighted in bold.

1) COMPARISON RESULTS FOR BAP
As shown in Table 8, the proposed BMO-HC performed better than 2 out of 5 algorithms and produced results same as the remaining algorithms in which they match the best-known solutions across all tested BAP instances. TS is the worse as it matched the best-known solutions on 13 out of 30 instances, while CG comes next as it matched 27 out of 30 instances. That means the BMO-HC results are very competitive and stable across all instances. This can be interpreted by the improved exploitation ability of BMO-HC which helps to search for high-performance regions of search space. Table 9 reports the comparison results of percentage deviation according to solution quality between BMO-SA and the following approaches: discrete spider monkey optimization (DSMO) [45], swap sequence based Artificial Bee Colony algorithm (ABCSS) [46], velocity tentative PSO (VTPSO) [47], GA [17], and PSO [17]. Table 9 demonstrate that the proposed BMO-SA algorithms performed better than GA and PSO algorithms across all tested instances, and better than the remaining approaches on 12 out of 15 tested instances. Even though DSMO perform better than other approaches on 3 instances, BMO-SA still in the first place where it outperformed all other approaches on 80% of the instances. In addition, BMO-SA matched the best-known solutions on 17 out of 20 tested TSP instances whilst none of the compared approaches has matched any of best-known solutions. This is due to the impact of SA algorithm in improving the search space and attaining wellbalance between exploration and exploitation of BMO.

IX. CONCLUSION
In this work, we proposed five hybrid approaches of BMO algorithm for solving combinatorial optimization problems: BMO-HC, BMO-LAHC, BMO-SA, BMO-IG, and BMO-VIG. The proposed approaches employed the single-based algorithms using low-level teamwork hybrid model.
Single-based algorithm used inside the BMO as a local search operator in which it searches the neighborhood of each solution in the population after each iteration. A hybridization probability is implemented to every solution before hybridizing to avoid early premature convergence. To evaluate the proposed methods, we tested them on the travelling salesman problem and berth allocation problem and all approaches are compared against each other. In addition, a comparative study was conducted against five recent methods in the literature including: GA, PSO, VTPSO, ABCSS, and DSMO. The computational results demonstrated that the hybrid BMO approaches were able to outperform the original BMO across all instances in both problem domains BAP and TSP. Furthermore, the BMO-HC and BMO-SA matched the best-known solutions.
across all instances for BAP and across most of TSP instances, respectively. Also, the hybrid BMO approaches showed very competitive results when compared to the stateof-the-art methods for BAP and TSP and ranked as among the best performing approaches. In conclusion, we proved that hybridizing BMO with single-based algorithms managed to improve its exploitation and attain the well-balance between exploitation and exploration of the algorithm.
In future work, it might be beneficial to investigate the possibility of proposing an online framework of the hybrid BMO to solve different combinatorial optimization problem without the need to select the single-based algorithm manually. On other words, proposing a BMO model that can automatically adjust the components of the hybrid BMO according to the search status.