An Improved Harris Hawks Optimization Algorithm Based on Bi-Goal Evolution and Multi-Leader Selection Strategy for Multi-Objective Optimization

2.1 Multi-Objective Optimization problem (MOP)

MOP involves optimizing problems using multiple objective functions [31]. This may be expressed as a minimization problem, as shown in the equation Eq. (1).

$\begin{aligned} & \operatorname{Min}: F_m(x), \\ & \left\{\begin{array}{l}g_i(x) \leq 0, i=1,2 \ldots, I \\ W_c(x)=0, c=1,2 \ldots, C \\ T_j \leq x_j \leq F_j, j=1,2 \ldots, n\end{array}\right.\end{aligned}$       (1)

The variable x represents a solution consisting of n decision variables, denoted by x₁, x₂, …, x_n subject to both inequality constraint "I " and equality constraints "C". "M" corresponds to the count of objective functions, while "T_j" and "F_j" represent the lower and upper bounds of the decision variable, respectively.

Definition (Pareto Dominance) Given two solutions "a" and "b", "a" is said Pareto dominates "b" (represented as a<b) if:

$\forall m \in\{1,2, \ldots, M\}: f_m(a) \leq f_m(b)$       (2)

$\exists m \in\{1,2, \ldots, M\}: f_m(a)<f_m(b)$       (3)

When no other feasible solution can dominate a particular solution, that solution is referred to as a non-dominated solution (Pareto solution). The collection of all non-dominated solutions is termed the Pareto set (PS), defined as:

$P S=\{a \in X \mid \nexists b \in X, a<b\}$       (4)

By applying the objective vector to the PS, we obtain what is referred to as the Pareto front (PF) representation.

2.2 Non-dominated sorting method (NDS)

The NDS [6] method, is a commonly utilized technique for solving multi-objective optimization problems. It involves selecting non-dominated individuals within the population and creating a hierarchical structure, representing a set of optimal solutions. The method sorts the population into several non-dominated layers, allowing for the identification of high-quality solutions that are diverse and well-spread across the pareto front. The algorithm assigns rank to individuals, which are determined by their degree of non-domination. The individuals that are best and not dominated by any other solution are given the first rank. This mechanism is repeated until the complete population has been sorted, and individuals are assigned ranks based on their level of domination. NDS is an effective method for finding a set of high-quality solutions that are well-spread across the pareto front. The principle of NDS is demonstrated in Figure 1.

1.png

Figure 1. Principle of Non-dominated sorting method [6]

2.3 Bi-Goal Evolution (BIGE)

Li et al. [30] proposed an innovative mechanism called Bi-Goal Evolution (BIGE). This mechanism is designed to formulate a selection strategy that is well-suited for extensive objective spaces. This approach converts a high-dimensional objective space into a bi-objective space, effectively considering the proximity and crowding degree of solutions within the population. The proximity calculation (Eprox) for a given individual s is computed by aggregating its individual objective values into a single comprehensive objective value, as shown in the Eq. (5):

$E_{\mathrm{Pr} o x}=\sum_{i=1}^k E^i(s)$       (5)

Here, “k” represents the total count of objectives, and Eⁱ(s) refers to the goal value of the individual "s" in the i^th goal.

Calculating the crowding degree (E_Crow) of an solution (s₁) within the pareto front set (F) involves the following estimation:

$E_{\text {Crow }}(s 1)=\left(\sum_{s 2 \in F ; s 2 \neq s 1} \operatorname{shar}(s 1, s 2)\right)^{1 / 2}$       (6)

where, shar(s₁, s₂) means sharing function between two individuals s₁ and s₂, which is elucidated in Eq. (7):

$\begin{aligned} & \operatorname{shar}\left(s_1, s_2\right)=\left\{\begin{array}{l}\left(0.5\left(1-\frac{\operatorname{Edist}\left(s_1, s_2\right)}{r}\right)\right)^2, \text { if } \operatorname{Edist}\left(s_1, s_2\right)<r, \\ E_{\operatorname{Pr} o x}\left(s_1\right)<E_{\operatorname{Pr} o x}\left(s_2\right) \\ \left(1.5\left(1-\frac{\operatorname{Edist}\left(s_1, s_2\right)}{r}\right)\right)^2, \text { if } \operatorname{Edist}\left(s_1, s_2\right)<r, \\ E_{\operatorname{Pr} o x}\left(s_1\right)>E_{\operatorname{Pr} o x}\left(s_2\right) \\ \text { rand }\left(), \text { if Edist }\left(s_1, s_2\right)<r, E_{\operatorname{Pr} o x}\left(s_1\right)=E_{\operatorname{Pr} o x}\left(s_2\right)\right. \\ 0, \text { otherwise }\end{array}\right.\end{aligned}$       (7)

In this context, Edist(s₁, s₂) is the euclidean distance between two solutions (s₁, s₂). The parameter "r" signifies the niche radius, computed using "N" and "K" which respectively denote the population size and the total count of objectives.

$r=\frac{1}{\sqrt[k]{N}}$       (8)

2.jpg

Figure 2. The principle of BIGE strategy [30]

More importantly, in this mechanism of BIGE, using the sharing function mentioned above permits neighbouring solutions to be placed at a distance. And then, the new parameters (E_Prox) and (E_Crow) are estimated using Eq. (5) and Eq. (6), respectively. The principle of BIGE is depicted in Figure 2.

2.4 Harris Hawks Optimization (HHO)

HHO [26] is a newly introduced Algorithm that emulates the Harris hawk chasing behaviour. Below, the mathematical model of this method is presented.

2.4.1 Exploration phase

During this stage, the Hawks are diffused in random areas and use two alternative strategies to find prey:

$\begin{aligned} & L(i t+1)=\left\{\begin{array}{l}L_{\text {rand }}(i t)-r_1\left|L_{\text {rand }}(i t)-2 r_2 L(i t)\right|, q \geq 0.5 \\ \left(L_{\mathrm{Pr} e \mathrm{y}}(i t)-L_{A v}(i t)\right)-r_3\left(T+r_4(F-T)\right), q<0.5\end{array}\right.\end{aligned}$       (9)

where, L(it+1) denotes the Location of Hawks for the subsequent repetition, L(it) is the actual location of hawks, L_Prey(it) is the actual location of the prey. r₁, r₂, r₃, r₄, and q are assigned randomly in (0, 1). F and T represent the upper and lower limits of positions, respectively, L_rand signifies an individual hawk selected at random from the population, and L_Av represents the mean locations of the individuals which is determined by:

$L_{A v}(i t)=\frac{1}{N} \sum_{i t=1}^N L_j(i t)$       (10)

In this equation, N, L_j  denote the total count of hawk individuals within the population, and the location of each hawk, respectively.

2.4.2 Transition phase

The HHO can exchange diversification and intensification stages in response to the fleeing energy of the rabbit (E), which is computed by:

$E=2 * E_0\left(1-\frac{i t}{i T}\right)$       (11)

where, E₀ refers to the rabbit’s energy at the starting point, which is choosing randomly from an interval of -1 to 1. If the value of E exceeds 1, the HHO prioritizes the diversification strategy. Conversely, if E is less than or equal to 1, the  intensification strategy is favored. iT is the maximum count of iterations.

2.4.3 Exploitation phase

Based on the prey’s fleeing energy (E) and their chance fleeing probability (r), the hawks employ four distinct pursuit ways. and switch between them as detailed below:

• Soft Besiege strategy {|E|≥0.5, and r≥0.5}

Hawks encircle the prey gently, causing it to become exhausted before launching a sudden attack, which is described as:

$L(i t+1)=\left(L_{\mathrm{Pr} e y}(i t)-L(i t)\right)-E \mid F^* L_{\mathrm{Pr} e y}(i t)-L(i t)$       (12)

The magnitude of the rabbit’s random jumps during the escape process can be expressed by F=(1-r₅)∗2, where r₅ represents a randomly generated value within the range of 0 to 1.

• Hard Besiege strategy {|E|<0.5, and r≥0.5}

When the rabbit reaches a state of exhaustion and lacks the energy to flee, it may become vulnerable to surprise attacks by hawks without needing them to encircle its intended target. This scenario is described as follows:

$L(i t+1)=L_{\mathrm{Pr} e y}(i t)-E \mid\left(L_{\mathrm{Pr} e y}(i t)-L(i t)\right)$       (13)

• Soft besiege with progressive rapid dives {|E|≥0.5, and r<0.5}

If the rabbit possesses sufficient energy to successfully flee, Hawks may employ a gentle encirclement strategy prior to launching a surprise attack. In such cases, the algorithm calculates the subsequent movement of the Hawks using the equation below:

$Y=L_{\mathrm{Pr} e y}(i t)-E\left|L_{\mathrm{Pr} e y}(i t)^* F-L(i t)\right|$       (14)

We presume that they will descend into a space of D dimensions, guided by the function of levy flight (LF) [32, 33], by adhering to the following formula:

$Z=Y+L F(D) * S$       (15)

We note that D pertains to the dimensions involved in the optimization issue, while S represents a randomly sized vector of 1*D. Furthermore, the calculation of the LF will be computed as follows:

$L F(l)=0.01 * \frac{u^* \sigma}{|v|^{(1 / \beta)}}$       (16)

u and v are selected at random from [0, 1], the value of β remains constant at 1.5, and σ is defined as below:

$\sigma=\left(\frac{\sin (\pi \beta / 2) * \Gamma(1+\beta)}{2^{\frac{(\beta-1)}{2} *} \beta^* \Gamma((1+\beta) / 2)}\right)^{\frac{1}{\beta}}$       (17)

Thus, the calculation model utilized for updating the Hawks’ positions is specified as follows:

$L(i t+1)= \begin{cases}\gamma, & \text { if } F(\gamma)<F(L(i t)) \\ \lambda, & \text { if } F(\lambda)<F(L(i t))\end{cases}$       (18)

• Hard besiege with progressive rapid dives {|E|<0.5, and r<0.5}

Hawks aim to minimize the gap between their mean locations and the fleeing rabbit. This which can be described as follows:

$L($ it +1$)= \begin{cases}\gamma, & \text { if } F\left(\gamma^{\prime}\right)<F(L(i t)) \\ \lambda, & \text { if } F\left(\lambda^{\prime}\right)<F(L(\text { it }))\end{cases}$       (19)

where,

$Y^{\prime}=L_{\mathrm{Pr} e y}(i t)-E\left|F^* L_{\mathrm{Pr} e y}(i t)-L_{A v}(i t)\right|$       (20)

$Z^{\prime}=Y^{\prime}+S^* L F(D)$       (21)

L_Av has been defined in Eq. (10).

This part evaluates the efficacy of the GMOHHO method in handling various MOPs. Firstly, an overview of the experimental setup is provided. Subsequently, the benchmark functions used and the parameters applied to the comparative methods in this study are described. Finally, the statistical findings of the evaluated algorithms are presented and analyzed.

4.1 Experimental setup

We evaluated and compared the GMOHHO algorithm with well-known algorithms from various perspectives using two performance indicators: the Inverted Generational Distance (IGD) [34] and the Hyper Volume (HV) [35]. GMOHHO were implemented using the MATLAB (R2016a) programming language on a 64-bit Windows 10 operating system. The experimental trials were executed using a computer that featured an Intel Core(TM) i5-7300HQ CPU and 8GB of RAM. To ensure reliable results and mitigate the impact of randomness, we performed 31 runs of each comparative method, with each run consisting of 1000 iterations. All test cases were conducted with a population size of 100.

4.2 Test functions and compared algorithms

Table 1. Attributes of test functions

Table 2. Parameters of algorithms

To evaluate the efficacy of the newly introduced GMOHHO algorithm, an extensive collection of benchmark functions was utilized for evaluation. These benchmarks were carefully selected from well-established categories, namely the ZDT-series [36] (consisting of five bi-objective functions), the DTLZ series [37] (consisting of a set of seven three-objective functions), and the MaF series [38] (encompassing seven many-objective functions). The attributes of these evaluation functions are outlined in Table 1.

For the purpose of conducting a comparative study, three reputable benchmarking techniques, namely MSSA [19], MOEAD [11], and MOGWO [14], were chosen. The initial parameters of each algorithm, as recommended in their original sources, are provided in Table 2. All test methods were conducted with a population volume of 100. to reduce the impact of randomness, we conducted 31 runs of each comparative method, each consisting of 1000 iterations.

4.3 Comparative experimental results on ZDT serie

In this section, we comprehensively evaluate the efficacy of the suggested GMOHHO algorithm specifically in the context of bi-objective problems represented by the ZDT series. To effectively assess its performance, we conduct a comparative analysis between the GMOHHO and other well-established multi-objective methods. The results obtained from this comparative analysis are presented in Tables 3 and 4, where we utilize the IGD and HV indicators, respectively, to measure the performance of the algorithm.

Table 3, which displays the IGD indicator, reveals that the GMOHHO algorithm consistently outperformed all other tested techniques across the entire range of ZDT benchmark problems. This outcome clearly signifies the superior performance of the GMOHHO algorithm in addressing the specific challenges posed by these problems. Similarly, Table 4, which presents the HV measure, consistently showcases the GMOHHO algorithm's superiority over the other methods across all test functions. The statistically acquired results not only validate the efficacy of the GMOHHO algorithm but also highlight its capability to effectively handle MOPs. However, to provide a more comprehensive analysis, it is crucial to delve into the reasons behind the observed performance differences among the algorithms and elucidate how the proposed GMOHHO method addresses any shortcomings of the other methods.

One potential reason for the GMOHHO algorithm's superior performance could be attributed to its innovative incorporation of hybridization between three strategies: the Bi-Goal Evolution (BIGE), Non-Dominated Sorting (NDS), and Multi-Leader Selection (MLS). This hybridization allows the GMOHHO algorithm to effectively balance exploration and exploitation, leading to improved convergence and search capabilities. In contrast, the comparative algorithms may exhibit limitations in terms of their exploration and exploitation capabilities, as well as their adaptability to diverse problem domains. These limitations can lead to suboptimal solutions and hinder their performance when faced with the complexities of MOPs.

It's essential to acknowledge that while the GMOHHO method demonstrated significant advantages, there are still potential areas for improvement. For instance, further investigations could be conducted to optimize the algorithm's parameter settings to enhance its overall performance. Additionally, the algorithm's robustness and scalability could be explored and validated on larger-scale MOPs to ensure its applicability to real-world problem domains.

Table 3. Results of IGD metric on ZDT test functions

The extensive evaluation of the GMOHHO in comparison to established methods on the ZDT benchmark series, as supported by both statistical metrics and visual convergence curves (Figures 5-7), confirms its superior performance in generating well-spread non-dominated solutions along the optimal pareto front.

Table 4. Results of HV metric on ZDT test functions

5.png

Figure 5. Best pareto fronts generated by each algorithm on ZDT1

6.png

Figure 6. Best pareto fronts generated by each method on ZDT2

7.png

Figure 7. Best pareto fronts generated by each method on ZDT4

4.4 Comparative experimental results on DTLZ serie

Table 5. Results of IGD metric on DTLZ test functions

Note: Best results are marked in bold

Tables 5-6 present the comprehensive statistical findings of the evaluated methods on the three-objective test problems from the DTLZ series, focusing on diversification evaluation criteria, convergence, and variety. The experimental findings demonstrate the superior performance of the proposed GMOHHO compared to other methods. This is evident in its enhanced convergence rate and its ability to generate a more diverse set of optimal solutions.

Table 6. Results of HV metric on DTLZ test functions

Note: Best results are marked in bold

8.png

Figure 8. Best pareto fronts generated by each method on DTLZ2

Table 5 summarizes the obtained findings of the IGD measure for the DTLZ series used in our study. The findings demonstrate that the GMOHHO method consistently achieves superior heterogeneity of non-dominated solutions compared to other well-known methods, as indicated by the lower IGD values. However, for the DTLZ2 and DTLZ6 test functions, the MOEAD method achieved slightly better results than our algorithm, ranking it second in these cases. Furthermore, Table 6 provides the findings obtained from the HV measure for the aforementioned benchmark tests. Our statistical analysis reveals that the optimal solutions generated by the GMOHHO outperform those of other well-known methods, particularly for DTLZ1, DTLZ2, DTLZ3, DTLZ5, and DTLZ7, as indicated by higher HV findings. However, it is noteworthy to mention that the MOGWO and MOEAD algorithms achieved better average HV values for DTLZ4 and DTLZ6, respectively.

Furthermore, we analyzed the convergence behavior of the the DTLZ series’ comparison algorithms and plotted the corresponding convergence curves, as shown in Figures 8-9. Our observations indicate that the optimal solutions obtained by the GMOHHO algorithm are well-distributed for the real pareto front for all the tested problems. Specifically, the suggested algorithm exhibited exemplary convergence behavior and achieved diverse Pareto solutions for all the DTLZ test problems. The algorithm's dynamic hybrid approach enables it to outperform other methods in most cases. However, certain areas for improvement have been identified, particularly regarding the performance on specific test functions. Future research can focus on addressing these limitations.

4.5 Comparative experimental results MaF on series

Table 7. Results of IGD metric on MaF test functions

Note: Best results are marked in bold

The results obtained from the benchmarks functions with many-objective within the MaF series, as presented in Tables 7-8, highlight the exceptional capabilities of the proposed GMOHHO in generating heightened convergence and a diverse set of optimal solutions. According to the IGD metric, the results in Table 7 demonstrate that the GMOHHO method consistently outperforms the competing methods across almost all the MaF series, converging quickly to the real pareto front. The only exception is the MaF1 test problem, where the MOEAD algorithm performed slightly better.

Table 8. Results of HV metric on MaF test functions

Analyzing the findings using the HV metric, as presented in Table 8, reveals that the GMOHHO algorithm surpasses the other methods for the majority of the examined issues, indicating its capacity to produce high-quality of solutions. However, similar to the IGD metric, the MOEAD algorithm exhibits a slight superiority over other algorithms for the MaF1 problem. Overall, these findings conclusively showcase the superior performance of the GMOHHO method in effectively addressing MOPs, surpassing the capabilities of the algorithms under investigation.

To gain deeper insights into the performance differences among the algorithms, let us discuss the underlying reasons. The GMOHHO algorithm employs a dynamic hybrid approach that combines BIGE, the NDS, and Multi-leader selection techniques. This unique combination empowers the algorithm to adeptly navigate the solution space, leading to a wider variety of optimal solutions. The method's efficient convergence to the real pareto front can be attributed to its ability to to harmonize diversification and intensification, thereby optimizing the trade-off between convergence and diversity.

Furthermore, Figures 10-11 provide graphical representations of the convergence behavior of the compared methods for the benchmark issues. The visualizations clearly illustrate that, across all the studied issues, the optimal solutions generated by the GMOHHO method are uniformly distributed and closely align with the ideal pareto front. This further substantiates the algorithm's ability to produce solutions characterized by a superior level of quality.

Despite the success of the GMOHHO method, it is important to acknowledge its limitations and identify areas for improvement. Further investigation could be conducted to understand the reasons behind the MOEAD algorithm's better performance on the MaF1 problem and explore potential modifications to enhance the GMOHHO algorithm's capabilities in solving such challenges. Additionally, future research could focus on evaluating our method on larger-scale and more intricacy multi-objective optimization issues to validate its scalability and robustness.

9.png

Figure 9. Best pareto fronts generated by each method on DTLZ7

10.png

Figure 10. Best pareto fronts generated by each algorithm on MaF1

11.png

Figure 11. Best pareto fronts generated by each method on MaF3

4.6 Running time analysis

The running time was recorded as the total computational time required by each algorithm to complete the optimization process on a given benchmark problem. This includes the time taken for initialization, iterations, and any necessary calculations. In our study, we conducted an evaluation of the algorithms' performance with respect to their running time on each test problem. To measure this, we recorded the time taken to complete an iteration of 1000. The outcomes are presented in Table 9.

Table 9. Execution time comparison of the four algorithms

Note: Best results are marked in bold

Upon analysis, we observed that the MSSA algorithm achieved the best overall results, ranking first for 12 test benchmarks. In terms of execution time, our proposed GMOHHO algorithm demonstrated the fastest performance for the DTLZ3, DTLZ4, DTLZ5, MaF3, MaF6, and MaF7 test problems. Our algorithm also achieved a second-place ranking for the remaining test functions, closely following the MSSA algorithm. There were a few exceptions, specifically for the ZDT2 and ZDT4 problems, where the GMOHHO algorithm ranked third behind the MOGWO method.

It is important to discuss the trade-offs between running time and efficiency indicators, such as IGD and HV. Although the GMOHHO algorithm might not consistently achieve the fastest execution time, it's essential to take into account the algorithm's overall effectiveness. When considering the quality of optimal solutions generated by the algorithms, as evaluated by IGD and HV metrics, the GMOHHO algorithm showcased competitive performance. It consistently outperformed the MOEAD and MOGWO algorithms. It is clear that there is a balance to be struck. Some algorithms may prioritize faster execution times at the expense of solution quality, while others may require more computational resources to attain improved performance. Regarding the GMOHHO, it strikes a favorable balance between running time and performance, achieving reasonable execution times while still generating high-quality non-dominated solutions.

4.7 Wilcoxon statistical test analysis

In this study, we incorporated the non-parametric Wilcoxon test [39] to assess the presence of significant statistical distinctions between the suggested GMOHHO algorithm and the other scrutinized algorithms being compared. The Wilcoxon test was chosen as a suitable statistical tool due to its ability to analyze paired data and handle non-normal distributions.

If the calculated p_value  is equal to or greater than the chosen significance level (α), the outcome would entail the non-rejection of the null hypothesis (H₀), indicating that there are no significant improvements associated with the comparative methods. However, if the p_valueis fewer than (α), we deny (H₀) to the benefit of the alternative hypothesis (H₁), which suggests that there are significant improvements associated with the comparative methods. In our analysis, we used a significance threshold of 5%, and compared the outcomes of two algorithms using the IGD and HV metrics. Specifically, µ1 and µ2 reflect the outcomes of the first and second algorithms, respectively.

The Wilcoxon analysis conducted on the IGD measure involved formulating the ensuing hypotheses:

• H₀: µ1≥µ2
• H₁: µ1<µ2

To perform the Wilcoxon analysis on the HV measure, we formulated the ensuing hypotheses:

• H₀: µ1≤µ2
• H₁: µ1>µ2

Table 10. Test Wilcoxon results for IGD metric

Note: Best results are marked in bold

The null hypothesis (H₀) assumes that the IGD and HV scores of the compared methods are either equal to or worse than those of the proposed GMOHHO. Conversely, the the opposing hypothesis (H₁) suggests that the IGD and HV scores of the GMOHHO algorithm are superior to those of the comparative methods.

Tables 10-11 present the outcoms of the Wilcoxon statistical evaluations, confirming that the alternative hypothesis (H₁) is accepted for the GMOHHO algorithm based on both the IGD and HV measures. The lower IGD values obtained by the GMOHHO algorithm indicate its superior ability to approach the true pareto front. Additionally, the higher HV values demonstrate its effectiveness in achieving better coverage of the Pareto set in contrast to the other methods, with a significance degree of less than 5%. These findings underscore the GMOHHO algorithm's potential in enhancing Pareto set diversity while converging towards the true Pareto front.

Table 11. Test Wilcoxon results for HV metric

Note: Best results are marked in bold

While the Wilcoxon test is a robust and widely used non-parametric test, it is important to acknowledge its limitations and assumptions. One assumption is that the observations within each pair are independent and identically distributed. Violation of this assumption may affect the interpretation of the results. Additionally, the Wilcoxon test assumes that the two samples being compared come from the same population, except for a shift in their distributions. It is essential to consider these limitations when interpreting the outcomes of the Wilcoxon test.

In summary, the non-parametric Wilcoxon test was employed to analyze the statistical significance of variances between the proposed GMOHHO and the comparative methods. Its selection was based on its suitability for paired data analysis and robustness against non-normal distributions.