A cluster-based co-evolutionary optimization method for bilevel multi-objective optimization

2.1. General framework of proposed algorithm

The framework of the CCBMO algorithm is established by algorithm 1. Initially, a set of $$N_u$$ upper level solutions $${\mathcal{X}}^u=\left\{{\mathbf{x}}_1^u,{\mathbf{x}}_2^u,...,{\mathbf{x}}_{N_u}^u\right\}$$ is randomly generated. Subsequently, the lower level optimization is conducted for $${\mathbf{x}}_i^u$$ individually via algorithm 2. The obtained lower level solutions and their corresponding upper level solutions are stored in $${\mathcal{L}}^p$$ (line 2). The solutions in $${\mathcal{L}}^p$$ are then ranked based on the upper level objective values and the degree of constraint violation, with the nondominated solutions being stored in $${\mathcal{X}}^o$$ (line 3).

Algorithm 1: Proposed algorithm

Algorithm 2: Lower level optimization

During the main loop, we execute the DE on $${\mathcal{X}}^u$$, with the resulting offspring stored in $${\mathcal{X}}^{uo}=\left\{{\mathbf{x}}_1^{uo},...,{\mathbf{x}}_{N_u}^{uo}\right\}$$ (line 5). For each solution in $${\mathcal{X}}^u$$, the closest solution to $${\mathbf{x}}_i^{uo}$$ is denoted as $${\mathbf{x}}_i^{um}$$ (lines 6-8). Next, we explain the several types of matches that can occur as a result. As shown in Figure 1, in the first case, $${\mathbf{x}}_1^{uo}$$ is matched with its closest parent upper level solution, $${\mathbf{x}}_1^u$$; therefore, $${\mathbf{x}}_1^{um}={\mathbf{x}}_1^u$$. In the second case, $${\mathbf{x}}_2^{uo}$$ and $${\mathbf{x}}_5^{uo}$$ are both matched with $${\mathbf{x}}_2^u$$, as they are the closest to it; therefore, $${\mathbf{x}}_2^{um}={\mathbf{x}}_2^u$$ and $${\mathbf{x}}_5^{um}={\mathbf{x}}_2^u$$; in the third case, no offspring are closest to the solutions of $${\mathbf{x}}_6^u$$ and $${\mathbf{x}}_7^u$$.

Figure 1. We provide an example of a block plot for matching and clustering in Figure 1A and B. (A) Clustering of $${\mathcal{X}}^{uo}$$ via kmeans; (B) match the solutions in $${\mathcal{X}}^u$$ with each solution in $${\mathcal{X}}^{uo}$$ based on the Euclidean distance.

Then $${\mathcal{X}}^{uo}$$ is clustered into $$N_u/2$$ groups using the k-means, denoted as $${\mathcal{G}}_1,\dots ,{\mathcal{G}}_{N_u/2}$$ (line 9). K-means is a clustering method based on Euclidean distance, which aims to minimize the sum of the distances of N objects from the nearest center point, and to divide these objects into K groups according to the size of the K value, and to consider the individuals in the same group to be similar. The index of each solution in $$i$$th group in $${\mathcal{X}}^{uo}$$ is stored in $${\mathcal{G}}_i$$. Following the execution of the k-means clustering, the solutions belonging to the same cluster exhibit proximity, suggesting similarities among the corresponding tasks of lower level optimization. Therefore, the lower level optimization is performed simultaneously on the solutions belonging to the same group based on algorithm 2. After that, the optimized lower level solutions and their corresponding upper level solutions are subsequently stored in $${\mathcal{L}}^s$$ (lines 11-14). The parents in $${\mathcal{L}}^p$$ and offspring in $${\mathcal{L}}^s$$ are combined, denoted as $$\mathcal{L}$$. Then, $$N_u$$ distinct upper level solutions are chosen to replace the original solutions in $${\mathcal{X}}^u$$, based on the nondominated sorting and crowding distance obtained from the upper level environmental selection. The selected upper level solutions and their corresponding lower level solutions are stored in $${\mathcal{L}}^p$$. Finally, the nondominated solutions in $$\mathcal{L}$$ are selected and merged into $${\mathcal{X}}^o$$ to update it.

In order to select the better solution for each iteration, the selection operator of NSGA-Ⅱ is applied to $$\mathcal{L}$$ to obtain $${\mathcal{X}}^u$$, and $${\mathcal{X}}^u$$ is used to store $$N_u$$ distinct upper level solutions. This selection is based on the upper level objective values and the upper level constraint violation. To be specific, the solutions in $${\mathcal{X}}^u$$ are initially divided into different nondominated sets using the constrained-domination principle. According to this principle, if we have two solutions, $${\mathbf{x}}_1$$ and $${\mathbf{x}}_2$$, $${\mathbf{x}}_1$$ is considered better than $${\mathbf{x}}_2$$ if any of the following conditions are met:

1. $$cv\left({\mathbf{x}}_1\right)=0$$ and $$cv\left({\mathbf{x}}_2\right)=0$$, and $${\mathbf{x}}_1$$ Pareto dominates $${\mathbf{x}}_2$$;

2. $$cv\left({\mathbf{x}}_1\right)=0$$ and $$cv\left({\mathbf{x}}_2\right)>0$$;

3. $$cv\left({\mathbf{x}}_1\right)>0$$ and $$cv\left({\mathbf{x}}_2\right)>0$$, and $$cv\left({\mathbf{x}}_1\right)<cv\left({\mathbf{x}}_2\right)$$.

Then, the best half of the upper level solutions are added to $${\mathcal{X}}^u$$.

2.2. Lower level optimization

algorithm 2 provides the process of lower level optimization. First, an empty set $$\mathcal{S}$$ is created (line 1). Then we use $$k$$ to denote the index of the $$j$$th element in $${\mathcal{G}}_i$$. Since we match a solution $${\mathbf{x}}_i^{um}$$ for each $${\mathbf{x}}_i^{uo}$$ (lines 6-8 of algorithm 1), we can find $${\mathbf{x}}_k^{um}$$ corresponding to $${\mathbf{x}}_k^{uo}$$. Then we store the lower level solution corresponding to $${\mathbf{x}}_k^{um}$$ into $$S$$. As shown in Figure 2, the $${\mathbf{x}}_1^{uo},...,{\mathbf{x}}_8^{uo}$$ were divided into four groups using the k-means algorithm. For example, $${\mathbf{x}}_1^{uo}$$ is in a separate group. Based on Figure 1, $${\mathbf{x}}_1^{uo}$$ is matched with $${\mathbf{x}}_1^u$$, and the lower level solutions whose upper level solution is $${\mathbf{x}}_1^u$$ are stored in $$S$$. $${\mathbf{x}}_2^{uo}$$ and $${\mathbf{x}}_3^{uo}$$ belong to the same group, but they are matched with solutions $${\mathbf{x}}_2^u$$ and $${\mathbf{x}}_3^u$$, respectively. Therefore, the lower level solutions whose upper level solutions are $${\mathbf{x}}_2^u$$ and $${\mathbf{x}}_3^u$$ are placed in $$S$$.

Figure 2. Comparison of the average FEs statistical results of the three algorithms on TP1-TP2 and DS1-DS5. FEs: Fitness evaluations.

After removing the repeat values in $$S$$, $$min\left\{N_l/2,\left|S\right|\right\}$$ solutions are randomly selected from $$S$$ and stored in $${\mathcal{X}}_1^l$$. Additionally, $$N_l-min\left\{N_l/2,\left|S\right|\right\}$$ lower level solutions are randomly generated and stored in $${\mathcal{X}}_2^l$$. Then, $${\mathcal{X}}_1^l$$ and $${\mathcal{X}}_2^l$$ are combined to generate $${\mathcal{X}}^l$$ (lines 7-10 of algorithm 2). During the main loop of algorithm 2, we execute the differential evolution on $${\mathcal{X}}^l$$, with the resulting offspring stored in $${\mathcal{X}}^{lo}$$ (line 12). We merge $${\mathcal{X}}^l$$ and $${\mathcal{X}}^{lo}$$, and store the combined set in $$\mathcal{C}$$. For solutions in $${\mathcal{G}}_i$$, we choose $$N_l$$ optimal lower level solutions from $$\mathcal{C}$$ considering the environmental selection of the upper level, and then store them in $${\mathcal{A}}_1,...,{\mathcal{A}}_{|{\mathcal{G}}_i|}$$. Eventually, the lower level solutions merge, replacing the original $${\mathcal{X}}^l$$ (lines 12-18). Upon meeting the termination conditions of the lower level optimization, we store $${\mathcal{A}}_1,...,{\mathcal{A}}_{|{\mathcal{G}}_i|}$$ and their corresponding upper level solutions in $${\mathcal{L}}_1^s,...,{\mathcal{L}}_{|{\mathcal{G}}_i|}^s$$.

2.3. Termination criterion

Termination conditions for both upper and lower layers are based on their hypervolume (HV). The rate of convergence of HV for both the upper and lower levels is called the H-metric and is calculated as follows:

where $$H{V}^{max}$$ and $$H{V}^{min}$$ denote the maximum and minimum values of the hypervolume in the output solutions, respectively, and $${}^{\prime }{u}^{\prime }$$, $${}^{\prime }{l}^{\prime }$$ are used to distinguish the upper and lower levels respectively. We use the maximum objective value of the nondominated solution as the reference point to calculate the HV.

When the number of iterations exceeds $$\sigma$$, and $$H<\epsilon$$, the algorithm terminates. We set the termination conditions as $${\epsilon }_l=0.001$$; $${\epsilon }_u=0.001$$; $${\sigma }_l=10$$, and $${\sigma }_u=40$$ dollars. In addition, for upper level optimization, if $${\sigma }_u$$ exceeds 40 and there are six consecutive generations where $$H_u$$ remains constant, the algorithm terminates. Since we optimize the lower level solutions simultaneously for multiple upper level solutions, for all lower level optimization processes, if $${\sigma }_l$$ exceeds 10 and there are five consecutive generations where the value of $$H_l$$ remains unchanged, the lower level optimization terminates.

3.1. Test problems

Two benchmark test sets, TP and DS ^[23], were selected for the experimental studies. In the TP test suite, we chose two test problems, TP1 and TP2, and the DS test suite includes five test problems (denoted as DS1-DS5).

The detailed settings for the TP test set and DS test set are provided in Table 1. $$D^u$$ and $$D^l$$ represent the dimensions of the upper and lower level decision variables, respectively. The dimensionality of both upper and lower level decision variables is K = 5 for DS1–DS3. For DS4-DS5, the dimensionality of the upper level decision variables is K = 1, and lower level decision variables is K + L = 10. The other parameters of the comparison algorithms remain the same as in their original papers.

3.2. Compared algorithms and parameter settings

In order to evaluate the performance of the algorithms, two algorithms were chosen for comparison: BLMOCC ^[20], MOBEA-DPL ^[27].

(1) BLMOCC: it is an algorithm that uses a knowledge-based variable decomposition strategy to solve a MBOP. The knowledge-based variable decomposition strategy is used throughout the optimization process, and the variables are divided into three groups according to the correlation between the variables in the two levels. Different optimization methods are used independently for different groups.

(2) MOBEA-DPL: it is an algorithm developed on the framework of the nested bilevel multi-objective optimization algorithm. It uses a dual-population optimization strategy to improve the solution of the lower level optimization problem. The first group is used to store nondominated solutions in the lower level, and the second group is used to store upper level solutions that are not dominated by solutions in the first group. Additionally, to increase the effectiveness of the search, the offspring of the upper level solutions are selected from the neighborhood of the best solutions in the existing solutions.

3.3. Performance metrics

We limit the maximum FEs of the upper and lower levels to facilitate comparing the number of true evaluations of different algorithms. The maximum number of evaluations for the upper and lower levels is 50, 000 and 1, 000, 000, respectively. $$F{E}^u$$ and $$F{E}^l$$ represent the number of evaluations in the upper and lower levels, respectively. In addition, we compare the sum of the actual number of evaluations in the upper and lower levels (i.e., $$F{E}^u+F{E}^l$$).

Without loss of generality, we consider using the inverted generational distance (IGD) ^[28] and hypervolume (HV) ^[29] to measure the performance of the algorithms. Both metrics can measure the diversity and convergence of the obtained solutions. A smaller IGD value means better algorithm performance, while the opposite is true for HV values, where a larger HV value means better algorithm performance.

3.4. Comparison between CCBMO and other algorithms on each index

In this case, we measure whether the algorithm can converge well with a small number of FEs. We recorded the average of FEs (including the upper level, lower level, and the sum of both levels), the IGD, and HV for all test problems. Specifically, each algorithm ran independently 21 times.

From Table 2, it can be seen that CCBMO has better results than BLMOCC and MOBEA-DPL in terms of lower level average FEs on TP1-TP2, DS1, and DS3-DS5, while slightly worse than BLMOCC in terms of upper level average FEs. This is because in lower level optimization, we introduce the idea of co-optimization, which optimizes similar lower level optimization problems at the same time, thus greatly reducing the number of optimizations and evaluations at the lower level level. However, as shown in Figure 3, due to the obvious advantage of CCBMO in the average FEs of the lower level, the sum of the average FEs of the upper and lower levels shows the best results, except that MOBEA-DPL outperformed CCBMO and BLMOCC on DS2. Obviously, this is all due to the improvement of the lower level optimization strategy.

Figure 3. Comparison of the statistical results of the average IGD values of the three algorithms on TP1-TP2 and DS1-DS5. IGD: Inverted generational distance

Tables 3 and 4 provide the comparative results of three algorithms on seven test problems in terms of IGD and HV. Figures 3 and 4 clearly demonstrate that CCBMO outperformed both BLMOCC and MOBEA-DPL. Next we analyze these experimental results in detail.

Figure 4. Comparison of the statistical results of the average HV value of the three algorithms on TP1-TP2 and DS1-DS5. HV: Hypervolume.

Both TP1 and TP2 are non-scalability problems. MOBEA-DPL outperformed other algorithms in terms of IGD value for TP1 due to its lower level dual-population strategy, but it inevitably increased the lower level FEs. CCBMO demonstrates a slightly inferior performance in terms of IGD value compared to MOBEA-DPL. Among all the compared algorithms, BLMOCC obtained the best HV value but the worst IGD value. Regarding TP2, CCBMO achieved the most satisfactory IGD and HV, followed by MOBEA-DPL. BLMOCC performed the worst.

DS1 is often used to evaluate algorithm performance and the ability of the algorithm to coordinate the processing of upper level and lower level tasks. Regarding DS1, MOBEA-DPL achieved the best IGD and HV. CCBMO performed slightly inferior to MOBEA-DPL in terms of IGD and HV values but significantly outperformed BLMOCC.

DS2 is used to assess the algorithm's capability to conduct extensive searches before converging on the frontier. For DS2, CCBMO significantly outperformed the contrasting algorithms. Compared with BLMOCC and MOBEA-DPL, CCBMO could reduce the IGD value by one to two orders of magnitude, respectively, and increase the HV value by two orders of magnitude. BLMOCC and MOBEA-DPL had difficulty achieving satisfactory results in the DS2 test problem.

DS3 involves discrete variables, and as the number of variables increases, the problem becomes increasingly challenging, making it difficult for traditional algorithms to address. MOBEA-DPL achieved the best IGD and HV on DS3. In contrast, CCBMO performed worse than MOBEA-DPL due to its utilization of a traditional nesting method for addressing this problem. BLMOCC had the worst IGD and HV.

Both DS4 and DS5 evaluate the algorithm's capability to search for the appropriate lower level frontier that corresponds to the upper level frontier. Additionally, DS4 necessitates identifying a specific point in the lower level solution that corresponds to the upper level solution. In the case of DS4, CCBMO achieved superior IGD and HV values, which are one and two orders of magnitude better than the corresponding IGD values of MOBEA-DPL and BLMOCC, respectively, and one order of magnitude better than the HV value of MOBEA-DPL. In the case of DS5, CCBMO had significant superiority over the comparison algorithms. In terms of IGD value, CCBMO outperformed both comparison algorithms by two orders of magnitude. The HV value of CCBMO is an order of magnitude better than that of MOBEA-DPL. Therefore, these experiments demonstrate the effectiveness of CCBMO in solving BMOPs.

3.5. Comparison of CCBMO with its variants

To illustrate the effectiveness of the mechanism in CCBMO, we compared two variants of CCBMO, CCBMO-M and CCBMO-C, where CCBMO-M represents the variant that removes the k-means clustering in CCBMO, and CCBMO-C represents the variant that removes the matching between the upper level offspring and parents in CCBMO. Tables 5-7 present the results of CCBMO, CCBMO-M and CCBMO-C regarding the average number of FEs, IGD and HV over 21 runs on TP1-TP2 and DS1-DS5, respectively.

Table 5 gives a comparison of CCBMO and its variants in terms of average FEs. CCBMO consumed the least number of evaluations in solving TP1-TP2, DS2 and DS4. CCBMO-M performed slightly worse than CCBMO, consuming the least number of evaluations on solving DS1 and DS3. CCBMO-C only consumed the least number of evaluations when solving DS5.

Tables 6 and 7 compare the IGD and HV values of CCBMO and its variants (CCBMO-M and CCBMO-C), respectively. In terms of TP1 and DS2-DS5, CCBMO outperformed CCBMO-M and CCBMO-C in both IGD and HV values. For TP2, CCBMO obtained the best IGD value, but the HV value was slightly worse than CCBMO-M. CCBMO-M had the best IGD and HV values when solving DS1. On the other hand, CCBMO-C performed the worst, not obtaining the best IGD and HV values for any problem. Therefore, our experiments demonstrate the effectiveness of matching and clustering incorporated in CCBMO.

Acknowledgments

Thanks to Yaochu Jin, Xingyi Zhang, Ran Cheng, Ye Tian, Xilu Wang, Dan Guo, Dimo Brockhoff, Handing Wang Qingfu Zhang and Chaoli Sun for providing their source code.

Authors' contributions

Proposal and design of research propositions, drafting or final revision of the paper: Hu WY, Huang PQ

Perform data collection: Hu WY

Technical and material support provided: Li F, Ge EQ

Availability of data and materials

Not applicable.

Financial support and sponsorship

The authors is supported by the National Natural Science Foundation of China under 61903003 and supported by Anhui Provincial Natural Science Foundation (Grant 2308085MF199 and 2008085QE227) and Scientific Research Projects in Colleges and Universities of Anhui Province (Grant 2023AH051124) and Supported by the Open Project of Anhui Province Engineering Laboratory of Intelligent Demolition Equipment (Grant APELIDE2022A007) and Supported by the Open Project of Anhui Province Key Laboratory of Special and Heavy Load Robot (Grant TZJQR001-2021).

Conflicts of interest

All authors declare that they are bound by confidentiality agreements that prevent them from disclosing their conflicts of interest in this work.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2024.