Cooperative multi-population Harris Hawks optimization for many-objective optimization

This paper presents an efficient cooperative multi-populations swarm intelligence algorithm based on the Harris Hawks optimization (HHO) algorithm, named CMPMO-HHO, to solve multi-/many-objective optimization problems. Specifically, this paper firstly proposes a novel cooperative multi-populations framework with dual elite selection named CMPMO/des. With four excellent strategies, namely the one-to-one correspondence framework between the optimization objectives and the subpopulations, the global archive for information exchange and cooperation among subpopulations, the logistic chaotic single-dimensional perturbation strategy, and the dual elite selection mechanism based on the fast non-dominated sorting and the reference point-based approach, CMPMO/des achieves considerably high performance on solutions convergence and diversity. Thereafter, in each subpopulation, HHO is used as the single objective optimizer for its impressive high performance. Notably, however, the proposed CMPMO/des framework can work with any other single objective optimizer without modification. We comprehensively evaluated the performance of CMPMO-HHO on 34 multi-objective and 19 many-objective benchmark problems and extensively compared it with 13 state-of-the-art multi/many-objective optimization algorithms, three variants of CMPMO-HHO, and a CMPMO/des based many-objective genetic algorithm named CMPMO-GA. The results show that by taking the advantages of the CMPMO/des framework, CMPMO-HHO achieves promising performance in solving multi/many-objective optimization problems.


Introduction
There are many practical multi-objective optimization problems (MOOPs) in engineering and scientific research [21,38,47]. A number of multi-objective evolutionary algorithms (MOEAs) and swarm intelligence optimization algorithms (MOSIOAs) have been proposed to solve MOOPs for they can obtain a well-distributed and well-converged set of near Pareto-optimal solutions. Although some traditional MOEAs [15,57] and MOSIOAs [11,41,51] can effectively solve twoand three-objective problems, they show poor performance Na Yang, Zhenzhou Tang and Xuebing Cai contributed equally overall to this work. B Qian Hu huqian@wzu.edu.cn 1 College of Computer Science and Artificial Intelligence, Wenzhou University, Wenzhou 325035, China 2 College of Computer and Software Engineering, Anhui institute of Information Technology, Anhui 241199, China in solutions diversity and convergence when solving manyobjective optimization problems (MaOPs). The reason is that their dominance rules are not sufficient to select the elite individuals for the next generation, nor to satisfy the requirement of the population size required for MaOPs [24,31]. Some MOEAs and MOSIOAs proposed later, such as NSGA-III [14], MOMPA [8], MOP-GSO [17] etc., can obtain wellconverged solutions, however, it is difficult to maintain solutions diversity when using only one population to optimize all objectives.
It has been proved that, compared with a single population, algorithms with multiple populations can effectively improve solution diversity for both single-objective optimization problems [5,32,33] and MOOPs [13,31,42] by dividing the population into multiple co-evolutionary/cooperative subpopulations. Nevertheless, some multi-population algorithms, such as MPEA/SG [30] and CIEMO/D [40], change the number of subpopulations during evolution, which makes it a tough challenge to determine the appropriate number of subpopulations when solving MaOPs. By contrast, in the multi-population multi-objective (MPMO) framework proposed in [50], the number of subpopulations is always equal to the number of optimization objectives. Each subpopulation only optimizes one objective and contains an equal number of individuals. However, the MPMO framework used the density-based selection operator to maintain population diversity, while ignoring that the computational expense of the density-based selection operator becomes considerably expensive with the increase of the number of objectives [14,25].
Motivated by the above, this paper presents a novel cooperative multi-population multi-objective swarm intelligent algorithm, known as CMPMO-HHO. Specifically, we first propose a novel cooperative multi-population multi/manyobjective framework with dual elite selection named CMP-MO/des. CMPMO/des has the features of independently optimizing one objective in one subpopulation, the cooperation among subpopulations in the global archive, the dual elite selection, and the logistic chaotic single-dimensional perturbation (LCSDP). These promising features help CMPMO-HHO dramatically promote the convergence and diversity of the obtained near Pareto-optimal solutions.
Secondly, we leverage the novel Harris hawks optimization (HHO) algorithm [1] as the optimizer of subpopulations. HHO is a state-of-the-art single objective optimizer that has the advantages of fast convergence, easy implementation, and friendly expansion, and has been applied to many practical applications. Although this work has chosen HHO as the optimizer, it should be noted that the proposed CMPMO/des framework can work with any other single-objective optimizer without modification.
We summarize the contributions of this paper as follows:

MOSIOAs and MOEAs with a single population
There have been many state-of-the-art MOSIOAs and MOEAs with a single population proposed to solve MOOPs. To name a few, NSGA-II [15] uses non-dominated sorting and distance-based method to select elite solutions. SPEA2 [57] leverages the environment selection strategy based on the nearest neighbors and achieves a good distribution on highdimensional optimization problems. CMOPSO [46] is a multi-objective particle swarm optimization (PSO) algorithm with fixed inertia weight. However, these algorithms show poor performance on MaOPs. Therefore, some single population algorithms were proposed to deal with MaOPs. In [14], NSGA-III was proposed which adopts the reference point (RP) based method to determine elite solutions. C-MOEA/D [2] is a decomposition-based evolutionary algorithm which extends the ability of MOEA/D to deal with constraints using adaptive constraint processing. In [27], the authors proposed a new bandit-based adaptive operator selection method over MOEA/D, namely FRRMAB, to automatically select appropriate operators in an online manner. hpaEA [7] first defined the nondominated solutions exhibiting evident tendencies toward the Pareto-optimal front as prominent solutions, using the hyperplane formed by their neighboring solutions, to further distinguish among nondominated solutions. Then, a novel environmental selection strategy was proposed to balance the convergence and diversity. In PREA [49], a strategy based on the parallel distance was introduced to select individuals in the promising region to ensure the population diversity. In [10], the authors put forward an enhanced version of NSGA-III, known as ANSGA-III, which takes the advantage of the adaptive RP. ar-MOEA [48] leveraged the preference angle and reference information-based dominance to solve MOOPs and MaOPs. RPD-NSGA-II [18] used a new decomposition-based dominance relation to deal with MaOPs and a new diversity factor based on the penaltybased boundary intersection method. In PICEA-g [43], the family of decision-maker preferences and candidate solution population co-evolve to optimize the targets.

MOSIOAs and MOEAs with multiple populations
The aforementioned works only involve one population. Considering the benefits introduced by multiple populations in solving MOOPs and MaOPs, multi-population algorithms have received more and more attention. A novel multipopulation evolutionary algorithm, known as MPMMOES, was proposed in [54] for solving multimodal MOOPs. The original population should be divided into two groups of subpopulations of equal size. One is designed to search for the optimal solutions in objective space. The other subpopulation focus to obtain high-quality optimal solutions in the decision space. The optimizer in [26] uses the dynamic population strategy. In [42], a random migration strategy in MOPSO was proposed to improve the diversity of the population. In [12], the Pareto envelope-based selection algorithm II (PESA-II) was proposed which uses the hyper-grids to keep the well-distributed solutions. Manzoor et al. [34] presented a multi-objective self-adaptive multi-population based Jaya algorithm (PMO-SAMP-Jaya) to optimally schedule the energy consumption in a smart building. A grid search based multi-population particle swarm optimization algorithm (GSMPSO-MM) was presented in [29] to handle multimodal MOOPs.
To take advantage of multiple populations more effectively, other researches have focused on the co-evolutionary optimization among them. Said et al. [37] proposed an Indicator-Based version of their recently proposed Co-Evolutionary Migration-Based Algorithm (CEMBA), named IB-CEMBA, to solve combinatorial multi-objective bi-Level optimization problems. In [20], the authors presented a new multi-population hybrid genetic algorithm (MPHGA) that combines the standard genetic algorithm with the alternative location and assignment algorithm. Wang et al. [44] proposed a dual-population based evolutionary algorithm. These two populations iteratively exchange the information obtained from elite solutions during the evolution to collaboratively search for the optimal solutions of the problem. Ben Mansour et al. [3] proposed a cooperative version of the multi-objective local search algorithm (IBMOLS) based on a quality indicator, called W -CMOLS.
The above-mentioned multi-population optimization algorithms were mainly designed for MOOPs. There have also been some works reported for MaOPs. Dai et al. [13] presented an improved evolutionary algorithm for solving the multi-objective optimization problems, which uses the improved K-dominance to rank the solutions in subpopulations to generate offsprings. Similarly, Zheng et al. [55] proposed an evolutionary algorithm based on M2M population decomposition and reference distance. For each subpopulation, the projection distance to the direction vector is used to enhance the selection pressure of Pareto dominance.
Co-evolutionary or cooperation among subpopulations has been also widely used in multi-population MOSIOAs and MOEAs for MaOPs. A coevolutionary particle swarm optimization with the bottleneck objective learning (BOL) strategy was proposed in [31] where all populations are coevolutionary in a distributed manner. Naidu [35] proposed a hybrid cooperative multi-objective invasive weed optimization (IWO) based on the space transformation search (STS). Chen et al. [6] proposed a novel multi-objective ant colony system based on the multi-objective co-evolutionary multi-population framework, in which two ant colonies are used to deal with the two objectives respectively. A novel interval multi-populations multi-objective optimization method called the interval cooperative multi-objective artificial bee colony algorithm (ICMOABC) based on interval credibility was proposed in [53] where the interval credibility is selected as the interval dominant method. Rakshit et al. [36] proposed a new MaOP algorithm that solves the MaOPs using an improved DE mutation strategy and optimizing objectives in parallel. Liu et al. [30] developed a multi-population evolutionary algorithm with a single-objective guide to tackle many-objective optimization problems. It exploits the merits of both multiple populations and single-objective optimization to balance diversity and convergence of the evolution process. Table 1 summarizes the aforementioned works.

Method
In this section, the details of the cooperative multi-population multi-objective Harris Hawks algorithm for many objectives are presented.

Double-evolved cooperative multi-population framework
To improve the diversity and convergence when solving MaOPs, the novel cooperative multi-population framework, named CMPMO/des, is proposed in this paper. The features of CMPMO/des are as follows (Fig. 1).
• Intra-subpopulation migration and optimization: In CM-PMO/des, each objective is randomly assigned to one subpopulation and all individuals are equally allocated to each subpopulation. Each subpopulation performs migration and optimization towards the only objective independently. This one-to-one correspondence framework facilitates its extension to MaOPs.
Obtain E i (t) by performing NDS and RP. 6: Obtain S hho i (t) by performing one iteration of the HHO upon E i (t). 7: Obtain S chaos (t) by performing LCSDP on S hho (t). 11:

12:
Obtain A(t) by performing NDS and RP on Q(t).

13: end for
The procedure of CMPMO/des is as follows: Step 1: Initialize the parent population P(0) with the size of M × N where M is the number of objectives and N is the size of each subpopulation. Create M subpopulations with the size of N , denoted as S 1 , S 2 , . . . , S M , where M is the number of optimization objectives. Create an empty archive A(0) = ∅.
Step 2: Step 3: In iteration t, perform one iteration of the HHO algorithm (see "HHO optimizer") upon E i (t). Then we obtain a new subpopulation with the size of N , denoted as S hho i (t), and a merged population with the size of M × N as Step 4 In iteration t, perform LCSDP (see "Perturbation strategies") on S hho (t) and we obtain another population with the size of M × N , denoted as S chaos (t).
Step 5: Select N elite individuals from Q(t) by performing NDS and RP and we obtain A(t) for the next iteration.
Step 6: If the termination condition is satisfied, CMP-MO/des stops; otherwise, goes to step Step 2.

HHO optimizer
Harris Hawks optimization (HHO) is a novel heuristic algorithm which solves optimization problems by simulating the rabbit-hunting behavior of Harris Hawks [1]. The process of HHO involves two phases, namely the exploration phase and the exploitation phase. In HHO, a Harris hawk represents a where E 0 is the initial energy generated randomly in (−1, 1), t is the current iteration, and T is the maximum number of iterations. HHO adopts different strategies according to E. Specifically, HHO is in exploration phase if |E| ≥ 1, otherwise, HHO is in the exploitation phase.
Exploration phase (|E| ≥ 1). In iteration t, each Harris hawk updates its position with two equal-opportunity strategies.
where X (t) is vector of the locations of all hawks in iteration t, x rd (t) is the location of a randomly selected Harris hawk, x rb (t) represents the location of the rabbit, r 1 , r 2 , r 3 , r 4 and q are random numbers in (0,1), which are updated in every iteration, ub and lb are the upper and lower bounds of the variable, respectively, and x a (t) is the average of all hawks position vectors for current population, which is calculated as follows: where X i (t) is the location of the ith hawk and N represents population size.
Exploitation phase (|E| < 1) In the exploitation phase, all hawks update positions with different strategies based on the random number r 5 in (0,1) and rabbits' escape energy E. If r 5 ≥ 0.5, |E| ≥ 0.5, Harris hawks take the action of soft besiege as follows: where X (t) = x rb (t) − X (t) is the difference between the hawks positions and that of rabbit in tth iteration, J = 2(1 − r 5 ) is the misleading jump strength in [0,2] when prey escaping. If r 5 ≥ 0.5, |E| < 0.5, Harris hawks take the action of hard besiege as follows: If r 5 < 0.5, |E| ≥ 0.5, all hawks adopt soft besiege with progressive rapid dives to update their locations as follows: where Y and Z are respectively calculated as follows: Algorithm 2 Procedure of one iteration in HHO Update position by (4). 7: else if r ≥ 0.5 and |E| < 0.5 then 8: Update position by (5). 9: else if r < 0.5 and |E| ≥ 0.5 then 10: Update position by (6). 11: else if r < 0.5 and |E| < 0.5 then 12: Update position by (11). 13: end if and where D is the dimension of the decision variables, S is a D-dimensional random vector and L F is the Lévy flight function based on the following rule: where u and v are random numbers in (0,1), β is default constant 1.5. If r 5 < 0.5, |E| < 0.5, Harris hawks perform hard besiege with progressive rapid dives to update their location as follows: where . Similar to soft besiege with progressive rapid dives, Y and Z are retained only when better fitness values are obtained.
The detailed procedure of one iteration in HHO is summarized in Algorithm 2.

Perturbation strategies
To further improve the local search ability and speed up the convergence of HHO, the logistic chaotic single-dimensional perturbation (LCSDP) strategy is introduced in this paper. Compared with the random search, LCSDP is more evenly Furthermore, it keeps dimension information of the optimal solutions by searching for solutions disturbed in a single dimension [28]. Thus, LCSDP has good pseudorandomness, ergodicity, and sensitivity to the initial values. Specifically, for each individual in S hho (t), the value of a randomly selected decision variable is perturbed as follows.
where i is the index of the randomly selected decision variable, ub i and lb i are the the upper and lower bounds of the ith decision variable, C(t) is the logistic chaotic sequence in iteration t which is iteratively updated as follows.
where μ is the logistic chaotic sequence parameter and the initial value C(1) is a random number in (0, 1) . The complete flowchart of CMPMO-HHO is shown in Fig. 2.

Algorithm complexity analysis of CMPMO-HHO in one iteration
We analyze the time complexity of CMPMO-HHO based on the CMPMO-HHO scheme presented in Algorithm 1. In summary, the overall CMPMO-HHO time complexity is O(4(2M N + N )N 2 ).

Experimental settings
All algorithms in this paper run on MATLAB R2020b with Intel (R) Core (TM) I7-9700 CPU@3.00 GHZ and Windows 10 operating system. The experimental results of all the comparison algorithms are obtained by PlatEMO [39]. The size of each subpopulation is set to 100.

Benchmark functions and comparative algorithms
To evaluate the performance of CMPMO-HHO on MOOPs with only two or three objectives, 34 popular benchmark functions were selected from five popular benchmark suites on MOOPs, including five ZDT [56] functions, nine WFG [22] functions, seven DTLZ [16] functions, ten UF [52] functions, and three LSMOP [9] functions. The maximal number of ZDT functions evaluations is 300, that of other benchmark function is 3000. Other details about the benchmark functions are presented in Table 2.
To evaluate the performance of CMPMO-HHO on MaOPs with 5 and 10 objectives, 19 benchmark functions, i.e., seven DTLZ functions, nine WFG functions, and three LSMOP functions are used. The maximal number of evaluations is 3000. The number of objectives and dimensions follows the recommendation in [9,22,56]. Table 3 shows the detail settings.
Moreover, to quantitatively evaluate the importance of the different strategies leveraged in the CMPMO/des framework, we separately evaluated the performances of the three variants of CMPMO-HHO, namely MOHHO/SP (multiobjective HHO algorithm with single population), CMPMO-HHO/NoP (CMPMO-HHO without LCSDP) and CMPMO-HHO/SES (CMPMO-HHO with single elite selection), and CMPMO-GA which is a CMPMO/des based many-objective GA.

Performance metrics
In this paper, the inverted generational distance (IGD) [4] is adopted as the performance metric for it can concurrently quantify the convergence and diversity of MOOPs algorithms. The smaller the IGD value is, the better the convergence and distribution of the algorithm are. Assuming that an algorithm obtains a set A of non-dominated solutions and gives a uniformly sampled reference points set P on real PF, the calculation of IGD(A, P) is as follows.
where |P| is the size of set P, . 2 represents the Euclidean distance in objective space. Each benchmark function is tested 30 times independently to avoid contingency, and the average value (AVG.) and variance (STD.) are taken as the final performance indicators. The best results are shown in bold.
Besides, the Wilcoxon signed-rank test [19] at 0.05 significance level is conducted in this paper to show the statistically significant differences between CMPMO-HHO and other popular MOEAs and MOSIOAs on each benchmark function. The true hypothesis of the Wilcoxon signed-rank test is the P value. There were statistically significant differences between the two algorithms if P < 0.05, otherwise there were no significant differences. Three symbols, namely +, −, and =, indicate the results of CMPMO-HHO is significantly better, worse, or equivalent to the corresponding competitor, respectively. More specifically, if P ≥ 0.05, symbol = is used, if P < 0.05, and IGD value of CMPMO-HHO is smaller than the algorithm used for comparison, it is indicated by the symbol +, otherwise by symbol −. Tables 4,5, and 6 present the IGD values of CMPMO-HHO and the comparative algorithms on MOOPs with two or three objectives and MaOPs with 5 and 10 objectives, respectively. We ranked them according to the IGD values obtained from each test function. The results show that CMPMO-HHO ranks first on 15 of 34 functions (including six WFG functions, two ZDT functions, two DTLZ functions, four UF functions, and LSMOP2) and ranks top two on 22 of 34 functions, which indicates that CMPMO-HHO is quite competitive on MOOPs. For MaOPs, CMPMO-HHO performs the best on 6 WFG functions and 2 DTLZ functions with 5 five objectives, and 5 WFG functions, 2 DTLZ functions, and LSMOP2 with 10 objectives. Notably, although there are not so many benchmark functions in which CMPMO-HHO ranks first, CMPMO-HHO ranks top two in 11 of all 19 test benchmark functions with five objectives, and 9 benchmark functions with ten objectives.

Experimental results and analysis
We can also observe that, benefiting from the novel CMPMO/des framework, CMPMO-GA significantly outperforms NSGA-II/III, which confirms the effectiveness of the CMPMO/des framework. On the other hand, although CMPMO-GA always ranks close to CMPMO-HHO, it is inferior to CMPMO-HHO in most cases. This is merely because HHO performs much better than GA due to its much more complicated migration mechanisms [1]. Moreover, the three variants of CMPMO-HHO are inferior to CMPMO-HHO, which implies that the multi-population mechanism, dual elite selection and LCSDP can effectively improve the performance of the CMPMO/des framework. Table 7 summarizes the overall scores and the rankings of CMPMO-HHO and the comparative algorithms in terms of IGD values. The overall score of each algorithm was derived from its rankings on all benchmark functions (as shown in Tables 4, 5 and 6). Concretely, the first ranked algorithm gets a M is the number of objectives and D is the dimension of decision variables one score, the second gets two scores, and so on. The overall score of each algorithm is the sum of its scores on all benchmark functions. Finally, we rank the algorithms according to the overall score. The smaller the overall score, the higher the ranking. It can be found that CMPMO-HHO always ranks first in solving both MOOPs and MaOPs, which verifies that       CMPMO-HHO has significantly high performance in optimizing multi-objective and many-objective problems. Thereafter, Wilcoxon signed-rank tests were performed to investigate whether CMPMO-HHO and the compared algorithms have statistically significant differences. The details are shown in Tables 8, 9, and 10 for MOOPs and MaOPs (five and ten objectives) experiments, respectively. Both three tables show that most of the P-values of the Wilcoxon signedrank test are less than 0.05, which demonstrates significant differences between CMPMO-HHO and other algorithms in a statistical sense. The last row of tables (highlighted in bold) summarizes the number of instances on each significance case (+/ = /−). It can be seen that the numbers of + symbols in both tables are always in the majority, indicating that CMPMO-HHO is always statistically superior to the other algorithms compared.
To visually demonstrate the differences between the PFs obtained by CMPMO-HHO and the real PFs, we draw comparisons of PFs in Fig. 3. Note that PFs of MaOPs are difficult to be visualized due to their high objectives and dimensions, therefore, only PFs on two and three objectives are presented in Fig. 3. The red diamonds represent the PF values obtained by CMPMO-HHO, and the black solid points represent the real PF which are provided by PlatEMO [39].
It can be observed the PFs obtained by CMPMO-HHO and the real PFs are highly coincident and the solutions found by CMPMO-HHO are uniformly distributed in the PFs in most of the benchmark functions. However, the PFs of UF4 and LSMOP2 obtained by CMPMO-HHO are well-distributed, but poorly converged. On the contrary, the PFs of DTLZ1, DTLZ7, UF1, UF8, UF9, and UF10 convergence well to the real ones, but their distributions are not uniform enough. In the DTLZ3, UF3, UF5, UF6, LSMOP1, and LSMOP3 functions, the performance of CMPMO-HHO is not satisfactory in both distribution and convergence. Moreover, it should be noted that the real PF of WFG3 provided by PlatEMO [39] does not cover the whole PF, since the PF of WFG3 has a nondegenerate part as well as the intended degenerate part [23]. Authors in [23] derived the PF of WFG3 which is well coincident with ours.
By analyzing the characteristics of these benchmark functions with relatively poor distribution or convergence, we can find that most of them are multimodal or with irregular PFs, which indicates that the CMPMO-HHO algorithm needs to be improved in dealing with such kind of problems in the future. Even so, we can observe that CMPMO-HHO not only ranks first overall in the tests on 34 MOOPs and 19 MaOPs (five and ten objectives, respectively), but also ranks first in benchmark problems that are multimodal or with irregular PFs, as summarized in Table 11, which strongly proves that CMPMO-HHO is considerably competitive. On the other hand, the fact that CMPMO-HHO does not achieve ideal PFs in all benchmark functions is also reasonable according to the no free lunch (NFL) theorem [45].      Fig. 3 Comparisons between the real PFs (provided by PlatEMO [39]) and PFs obtained by CMPMO-HHO. † Note that the real PF of WFG3 does not cover the whole PF, since the PF of WFG3 has a nondegenerate part as well as the intended degenerate part [23]

Conclusion
This paper presents a swarm intelligent algorithm based on HHO, named CMPMO-HHO, to deal with multi-objective and many-objective problems. To ensure scalability over optimization objectives, CMPMO-HHO uses a novel cooperative multi-population framework named CMPMO/des. CMPMO/des includes four excellent strategies, namely the one-to-one correspondence framework between the optimization objectives and the subpopulations, the global archive for information exchange and cooperation among subpopulations, the LCSDP strategy, and the dual elite selec-tion mechanism based on the fast non-dominated sorting and the reference point-based approach, which help CMPMO/des achieve considerably high performance on solutions convergence and diversity.
We have conducted extensive comparative experiments on 34 multi-objective and 19 many-objective popular benchmark problems to evaluate the performance of CMPMO-HHO. The compared algorithms include 13 state-of-the-art MOEAs/MOSIOAs, three variants of CMPMO-HHO and a CMPMO/des based many-objective GA. The experimental results show that by taking the advantages of the CMPMO/des framework, CMPMO-HHO achieves an amaz- Table 11 Rankings of all algorithms in multi-objective test problems that are multimodal or with irregular PFs