Traffic flow prediction with a multi-dimensional feature input: A new method based on attention mechanisms

1.   Introduction

Multi-objective optimization problems (MOPs) exist in many fields of production practice ^[1]. Thus, lots of multi-objective evolutionary algorithms (MOEAs) have been designed to deal with MOPs with two or three objectives ^[2]. However, with the continuous improvement of people's requirements for sophistication, optimization models established in many complex industrial processes ^[3,4] and industrial applications ^[5,6,7] often involve more than three objectives, which are called as many-objective optimization problems (MaOPs). Generally, an MaOP involving $m$ conflicting objectives can be formulated as follows:

where $x = {({x_1}, {x_2}, \ldots, {x_n})^T}$ is the decision vector in the decision space $\Omega$. $F(x)$ denotes the objective vector and $f_i(x)$ is the $i$th objective ($i = 1, 2, \ldots, m$). Due to the conflicting property of different objectives, it is almost impossible to find a solution that optimizes all objectives simultaneously ^[8]. On the contrary, a set of trade-off solutions is needed, which is known as the Pareto front (PF) in the objective space and the Pareto set (PS) in the decision space ^[9].

Since the number of optimization objectives are increased, most of the solutions based on Pareto-dominance method become nondominated ^[10] and considerably slows down the evolutionary process, which leads to most Pareto dominance based MOEAs that fail to solve MaOPs. Besides, as the dimension of the objective space increases, the finite solutions become increasingly scarce in the high-dimensional objective space. As a result, many frequently used diversity maintenance strategy, such as crowding distance ^[11], are no longer effective on MaOPs.

To better solve MaOPs, researchers have proposed many different methods, which can be divided into the following five categories ^[12]:

The first category is modifying Pareto-dominance relation or developing new dominance relations to increase the selection pressure toward the PF for MaOPs. Examples of modifying dominance definitions include $k$-optimality relation ^[13], average and maximum ranking relations ^[14], $\epsilon$-dominance ^[15], $L$-optimality ^[16], grid-dominance ^[17], preference order ranking ^[18], fuzzy dominance ^[19], $\theta$-dominance ^[20] and $RP$-dominance ^[21]. The distinctive modifications have shown to be very promising for MaOPs. However, most of them can only be used within a specific algorithmic framework. Another typical idea in this category is to adjust the selection pressure through changing the dominated area of each solution, such as $\alpha$-dominance ^[22]. Two representatives are $CE$-dominance ^[23] and $CN$-dominance ^[24]. The former can dynamically adjuste the solutions' dominant area by varying a parameter $S$, and thereby adjust the selection pressure. The latter transforms the original objectives by using the author designed nonlinear function to further adjust the area dominated by each solution. While the idea of dynamically adjusting the dominated area of each solution enables Pareto-based MOEAs to enhance selection pressure among solutions, ensuring the quality of the selected solutions remains a challenging endeavor.

The second category is the indicator-based approach, which uses an evaluation indicator instead of Pareto-dominance to evaluate the quality of solutions ^[25]. Among the current indicators available, the IGD indicator as a comprehensive indicator is used to select the potential solutions in each generation (MaOEA-IGD) ^[26]. The hypervolume (HV) indicator ^[27] possesses good theoretical properties and is often used as the indicator function in multiobjective search. However, the high computation cost of the hypervolume holds back its spread, especially for MaOPs with more than five objectives ^[28]. Bader and Zitzler ^[29] suggested a new HV-based algorithm (HypE), which adopts Monte Carlo simulation to replace the exact hypervolume calculations, causing the computational complexity been alleviated, but in turn, the inexactly calculation of HV deteriorates the performance of HypE on MaOPs ^[30].

The third category is the objective reduction-based methods. Since some many-objective optimization problems have redundancy between the objectives, this method aims to identify the most relevant objectives first and then expurgate redundant ones that do not contribute significantly to describing the Pareto front. Deb and Saxena ^[31] developed a novel algorithm called PCA-NSGA-II, which combines principal component analysis (PCA) with NSGA-II for effectively handling MaOPs that have redundant objectives. Moreover, Singh et al. ^[32] designed a novel approach that emphasizes a reduced set of objectives instead of dealing with the complete dimensionality of the MaOP. Meanwhile, they adopt corner solutions, which are located in the boundaries of PF, to estimate the dimensionality of the true PF. Recent research has confirmed that certain methods based on objective reduction are susceptible when dealing with optimization problems that involve a high-dimensional PF ^[33]. In other words, these techniques may be effective for addressing problems that involve a moderate number of conflicting objectives.

The fourth technique is to use preference incorporation-based approaches, where the timing of integrating preference information into the optimization process is a crucial factor. The preference information given by decision maker (DM) can make the search direction that is biased towards the area of the PF, which is focused by DM, the many-objective optimization algorithms combined with the preference information can reduce the consumption of computing resources during the search ^[33]. For example, Wang et al. ^[34] developed a novel algorithm named PICEA-g, which enables decision-makers to obtain a comprehensive and diverse representation of the Pareto fronts, prior to eliciting and applying their preferences. In ^[35], a reference vector guided evolutionary algorithm is proposed, where the reference vectors not only decompose a MaOP into single-objective subproblems, but also elucidate user preference to target a preferred subset of the whole PF. This shows promising performance on certain problem types, but they suffer from suboptimal results since the bias may lead to premature convergence ^[36].

The last technique is a decomposition-based method, which decomposes a complex MOP into a number of subproblems and simultaneously optimizes them in a collaborative manner, is often known as another promising way. MOEA/D ^[37] is the most typical representative of this method. By assigning a set of evenly distributed weight vectors, MOEA/D maintains population diversity better and has the ability to solve many kinds of optimization problems with varying degrees of success ^{[38,39,40,41,42,43,44]}. Based on the framework of MOEA/D, MaOACO-RP utilizes the designed penalty boundary intersection and adaptive reference points to pick out solutions for the next generation ^[38]. Zhang et al. ^[45] integrated the designed information feedback model into MOEA/D for solving MaOPs with many decision variables. MaOEA-IT ^[46] adopts two stages to balance the convergence and diversity. Recently, the reference point-based MOEAs are highlighted in many-objective optimizations. A little bit different from decomposition-based MOEAs, the reference points are used to facilitate the diversity of selected solutions by associating each solution in objective space with them, e.g., NSGA-III ^[47], which adopts a reference point-based niche-preservation operation to replace the crowding distance operator in NSGA-II. Recently, Yi et al. ^[48] combines the crossover operator with the NSGA-III for MaOPs with many decision variables. However, with the increasing number of objectives, Pareto-dominance gradually loses its ability to filter solutions, thus, the solutions obtained by NSGA-III stress diversity more than convergence. In literature ^[20], a new algorithm called $\theta$-DEA is proposed to solve MaOPs, where a reference points-based dominance relation ($\theta$-dominance) is designed to select the solutions. The experimental results show that $\theta$-DEA can better improve the performance of NSGA-III for MaOPs, especially the convergence. In SPEA/R ^[30], it designs a new reference direction based fitness assignment scheme and an environmental selection strategy to improve the performance of SPEA2 in solving both multi-objective and many-objective problems, but the performance of SPEA/R remains constrained to Pareto-dominance with the number of objectives increasing.

In order to better maintain the diversity of solutions in objective space, and motivated by the idea of dividing objective space by well-distributed reference vectors to select the quality of obtained solutions, we propose a double association-based evolutionary algorithm for solving MaOPs. In comparison with existing approaches that rely on reference points/vectors, this paper offers the following key contributions:

$1)$ A double association strategy is designed to associate each solution with a reference vector. When the objective space is divided by the well-distributed reference vectors, different from other existing association methods ^[20,30], the double association strategy takes the empty subspace into account and associates it with the nearest solution. This new association method can facilitate the exploration of unknown areas.

$2)$ A new quality evaluation scheme is proposed to quantify the quality of each solution in subspace, where the convergence and diversity of each solution is first computed, and then a dynamic penalty coefficient is designed to balance the convergence and diversity by penalizing the diversity distribution of solutions.

The remaining structure of this paper is as follows: Section 2 elaborates the 3DEA in detail. Section 3 presents the experimental design for solving many-objective optimization. In Section 4, the experimental results are conducted and analyzed. Finally, Section 5 presents the conclusions of the study.

2.   Proposed algorithm: DAEA

2.1.   General framework

The framework of DAEA is given in Algorithm 1 and Figure 1 shows the flowchart of DAEA, which mainly corresponds to Algorithm 1. First, we employ the Das and Dennis's ^[49] systematic approach to generate evenly reference vectors. Then, an initial population $P_{0}$ with $N$ individuals are randomly produced. The ideal point $z^{*}$ and nadir point $z^{nad}$ are initialized in step 3 and step 4, respectively. They are computed by the minimum value and maximum value of obtained population $P_{0}$ for objective $f_{i}$, respectively, and both of them are updated in each iteration. Steps 6–23 denote that DAEA enters the iterative procedure. In step 7, on the basis of the preserved parent population $P_{t}$, the offspring population $Q_{t}$ is produced by employing the widely used genetic operators, i.e., the simulated binary crossover (SBX) ^[50] and the polynomial mutation (PM) ^[51]. Meanwhile, a combined population $R_{t}$ is formed by combining the parent population $P_{t}$ and offspring population $Q_{t}$. The population $S_{t} = \bigcup^{\tau}_{i = 1}F_{i}$, where $F_{i}$ is the $i$th Pareto nondominated level of $R_{t}$ and $\tau$ satisfies $\sum^{\tau-1}_{i = 1}|F_{i}| < N$ and $\sum^{\tau}_{i = 1}|F_{i}| \geq N$. In step 11, the normalization procedure ^[20] is implemented to normalize the objective population $\{F(x)|x\in S_t\}$ of $S_{t}$ to $\tilde F(S_t) = \{\tilde F(x)|x\in S_t\}$, and then each member in $\tilde F(S_t)$ is associated with a reference vector by our proposed double association strategy in step 12. For each solution in its associated subspace, the proposed quality evaluation scheme is used to evaluate the quality of each solution, and then sort them according the evaluated values, {the first-ranked solution in each subspace are put into the set $T(1)$, the second-ranked solution in each subspace are put into the set $T(2)$, and so on. Then steps 14–21 fill the population slots in $P_{t+1}$ using one level at a time, starting from $T(1)$, for the solutions in last accepted level $T(i)$, we randomly select the number of needed solutions in $T(i)$, since our proposed quality evaluation scheme takes both convergence and diversity into account. In the following sections, the implementation details of important component of DAEA will be described step by step.

2.2.   Reference vectors generation

In order to maintain the diversity of obtained solutions, evenly reference vectors $W = ({\lambda _1}, {\lambda _2}, \cdots, {\lambda _N})$ are generated by Das and Dennis's ^[42] systematic approach in DAEA, which divide the objective space into $N$ independent subspaces and reserve the most promising solutions of each subspace as much as possible. In this approach, reference vectors are sampled from a unit simplex and the number of reference vectors is equal to $N = \left(\begin{array}{c} H + m{\rm{ - 1}}\\ {\rm{ }}m{\rm{ - 1}} \end{array}\right)$, where $m$ is the dimension of objective space and $H$ is the number of divisions considered along each objective axis.

Suppose that ${\lambda _i} = (\lambda _i^1, \lambda _i^2, \cdots, \lambda _i^m)$ is the $i$th reference vector, let us consider

Then, the elements of reference vector ${\lambda _i}$ is obtained by

It can be clearly observed that when $H \geq m$, the above method will result in lots of reference vectors. Even if $H = m$, $N = \left(\begin{array}{c} 8 + 8{\rm{ - 1}}\\ {\rm{ }}8{\rm{ - 1}} \end{array} \right) = {\rm{6435}}$. This will significantly increase the computational cost. However, if we adopt the method reducing $H$ to relieve the computational burden, the generated reference vectors may unevenly distribute. To avoid such situation, the two-layered reference vectors with small values of $H$ as suggested in ^[47] is adopted in DAEA. Supposing the number of divisions of boundary and inner layers is $H_{1}$ and $H_{2}$, respectively, then the reference vectors are generated by Das and Dennis's method. The number of reference vectors (i.e., population size) is computed as

2.3.   Double association strategy

After the whole objective space is divided into $N$ independent subspaces, the population distributed in the target space needs to be assigned to different subspaces. The idea of assigning the population by associating each member in the normalized objective set with its closest reference vector has been employed in many recent papers, such as ^{[20,21,30,43,52]}, etc. However, they have different characteristics and motivations. In ^[52], when a subspace contains no solution, i.e., the subspace is empty, a solution is randomly chosen and associated with this subspace. This means a randomly selected solution is assigned to the empty subspace. This kind of association methods is the random association. Random association aims to increase the probability of the unknown region being explored, however, the operation of randomly selecting solutions may make the search even worse. In ^[34], each nondominated solution is associated with a subspace and then the density of each subspace is estimated by counting the number of solutions associated with the subspace. Since there may be some empty subspaces, it may lead to some incorrect estimation. Also, in ^[20,21,30], when empty subspace appears, these association methods have no further action and directly ignore them, which will destroy the diversity of the obtained solutions. This kind of association methods is the single association.

To better explore the unknown region and increase the diversity, we design a double association strategy, in which we do the first step of association as same as that in the existing one ^[53], i.e., each reference vector will generate a subspace, and each solution will be associated with its closest subspace, but like the discussion above, this can result in empty subspace(s) and can not keep the good diversity of the obtained solutions. To overcome this shortcoming, we do the second step of association. For each empty subspace, we associate a closest solution with it, in which the perpendicular Euclidean distance $d_{2}$ (the distance of each solution from each reference vector) ^[46] is used to measure the distance between the solution and the subspace (represented by the corresponding reference vector). In this way, the diversity of the obtained solutions will be greatly improved. The detail is as follows.

Currently, we have a normalized union population $\tilde F(R_t)$ (see Algorithm 1) in the objective space. Suppose that $\tilde F(x) = {({\tilde f_1}(x), ..., {\tilde f_m}(x))^T}\in \tilde F(R_t)$, the ideal point is the origin, then the Euclidean distance between the origin and the foot of the normal drawn from the solution $x$ to the reference vector ${\lambda _i}$ is denoted as $d_1^i(x)$, and the length of the normal is denoted as $d_2^i(x)$. To be specific, their mathematical descriptions are ^[16]:

For visually presenting them, Figure 2 illustrates the two distance measures $d_1^i(x)$ and $d_2^i(x)$ with respect to the reference vector ${\lambda _i}$ in a two-objective minimization problem.

In the first step association, for each solution $x_i\in R_t$, we first compute its perpendicular Euclidean distance $d_2^j({x_i})$ between $\tilde F({x_i})$ and each ${\lambda _j}(j = 1, 2, ..., N)$, then the solution $x_i$ is associated with its closest subspace (corresponding to the closest reference vector) based on the distance $d_2^j(x_i)$ for $j = 1, 2, ..., N$. After all solutions in $R_t$ have been associated with their closest subspaces, we begin the second step association. We examine all subspaces to see whether there is an empty subspace. If yes, we associate the empty subspace with its closest solution in $R_t$ using perpendicular Euclidean distance. Otherwise, we do not use the second step association. In this way, not only can the diversity of the obtained solutions be improved, but also the solution which is closest to the empty subspace will be preserved, which is very helpful to improve the diversity and can generate potential solutions in the subsequent evolution. The detailed description of the double association strategy is given in Algorithm 2.

To intuitively demonstrate the advantages of our proposed double association strategy, we compare it with two kinds of existing association schemes: Random association (randomly chooses a solution to associate it with an empty subspace in the second step association), e.g., ^[52], and single association (without second step association), e.g., ^[20,21,30]. Figure 3 gives a simple example to illustrate the influence of the random association, the single association and the proposed double association strategy on the performance of the algorithm. For these three association methods, the first step association operation of them are the same, we suppose that after the first step association operation, solution $a$ is associated with sub-space $S(1)$ and $b, c, d$ are all associated with sub-space $S(2)$. There is no solution associated with sub-space $S(3)$. If one uses the single association, such as ones in ^[20,21,30], $S(3)$ will be empty. This will result in bad diversity of solutions. If one uses the random association, a randomly chosen solution among $\{a, b, c, d\}$ (such as $d$) is assigned to $S(3)$, which also results in bad diversity of the solutions. If one uses our double association, since $c$ is closest to ${\lambda_3}$, it will be directly assigned to associate with subspace $S(3)$, which helps to improve the diversity of the solutions. This is very necessary for exploring unknown areas. Besides, the survival of solution $a$ in subspace $S(1)$ is very important for the unexploited area.

2.4.   Quality evaluation scheme

For comprehensively assessing the quality of each solution, we design a novel quality evaluation function for the solutions in each subspace $S(i)$, where this evaluation function takes into account both the convergence measure and diversity measure of each solution. The mathematical description of the new quality value of solution $x$ is as follows:

where $d_{1}^{i}(x)$ means the shortest perpendicular Euclidean distance of solution $x$ from the ideal point and is used to measure the convergence of solution $x$. The smaller the value of $d_{1}^{i}(x)$, the better convergence of solution $x$. The shortest perpendicular Euclidean distance $d_{2}^{i}(x)$ among the distances of each solution $x$ from each reference vector $\lambda_i$ is used to measure the contribution of solution $x$ to the global diversity, the smaller the value of $d_{2}^{i}(x)$, the better the diversity of solution $x$. $\mathop {\min }\limits_{y \in S(i), x\neq y} d(x, y)$ represents the minimum distance from the solution $x$ to all other solutions associated to the subspace $S(i)$, which is used to measure the contribution of solution $x$ to the local diversity. For diversity measure, the smaller the value of $d_{2}^{i}(x)$ and $\mathop {\min }\limits_{y \in S(i), x\neq y} d(x, y)$, the better the diversity of solution $x$. The parameter $\mu$ is a penalty parameter, and is designed to punish the solutions in dense area from global perspective. Considering that in an ideal state, each subspace should have two solutions distributed, here we set $\mu = \frac{|S(i)|}{2}$. As the subspace $S(i)$ is distributed with more solutions, the parameter $\mu$ is going to get bigger, which means stronger punishment will be inflicted on local diversity. Therefore, when a solution obtains a smaller quality value $FV$, it indicates that the solution possesses superior convergence and diversity characteristics.

Therefore, the new quality evaluation simultanously consider convergence measure and diversity measure and can adaptively protect the solutions located in the sparse areas. The quality value of a solution is lower, which denotes the solution has better convergence and diversity. After the solutions in each subspace are assigned to the corresponding quality value, we ascending sort the solutions in subspace $S(i)$ as $Num_{i}{(1)}, ..., Num_{i}{(|S(i)|)}$ based on the quality value $FV$. Then, the first-ranked solution in each subspace are put into the set $T(1)$, the second-ranked solution in each subspace are put into the set $T(2)$, and so on. The details are shown in Algorithm 3.

2.5.   Computational complexity of DAEA

The computational complexity of DAEA in one generation is dominated by the double association operator that is described in Algorithm 2. For the procedure of double association, there are two steps: In the first step of association (lines 2–8 of Algorithm 2), for a solution $x$, calculating $N$ perpendicular distances from $x$ to $N$ reference vectors needs $O(mN)$ computational cost, and calculating the shortest perpendicular distance from $x$ to $N$ reference vectors requires $O(N)$ computational cost. Thus, the association between $2N$ solutions with $N$ reference vectors needs $O(mN^2)$ computational cost. In the second step association (lines 9–14 of Algorithm 2), finding the empty subspace from $N$ subspaces and then associating it with a solution closest to it needs $O((N-1)N)$ operations in the worst situation. Thus, the overall worst complexity of one generation of DAEA is approximately $O(mN^2)$.

3.   Experimental design

In this section, experimental design is given to verify the performance of the proposed DAEA on many-objective optimization problems. First, the test problems and the compared algorithms are presented. Then, we will give a brief rundown of the performance metrics utilized in our experimental study. Next, we give parameter settings in this work. {All experiments are carried out in the open-source software PlatEMO ^[54] that has been widely used in many-objective optimization.

3.1.   Test problems

In order to test the performance of six algorithms involved in this paper, DTLZ test suite (DTLZ1-DTLZ7) ^[55] and WFG test suite (WFG1-WFG9) ^[56] as two common used many-objective optimization benchmarks are adopted in our experiments. They contains many different features, such as degenerate, biased, large scale, non-separable and partially separable of the decision variable in decision space as well as linear, convex, concave, mixed geometric structures and multi-modal of the PFs in objective space. These different features pose great challenges to the comprehensive performance of an algorithm. Besides, the objectives and decision variables of all these test problems can be scaled to any number.

In our experiment, the number of objectives for all test problems are taken from 5 to 20, i.e., $m \in \{5, 8, 12, 16, 20\}$. For DTLZ test suite, the number of decision variables is given by $n = m+k-1$. Different test problems may have different values of $k$. As suggested in ^[41,48], the parameter $k$ is set to 5 for DTLZ1, 10 for DTLZ2-DTLZ6 and 20 for DTLZ7. For WFG test suite, as suggested by the references ^[50] and ^[49], the number of decision variables for all test problems are set to 24 and the position related parameter is set to $m-1$.

To test the performance of our proposed DAEA in solving MaOPs, five well-known MaOEAs are selected to compared with DAEA based on a series of experiments. To be specific, NSGA-III ^[47] is selected since it is typically used as a baseline for solving MaOPs. $\theta$-DEA ^[20] adopts the single association method to decompose the population. SPEA/R ^[30] is mixed algorithm that combines decomposition and Pareto dominance. MaOEA-IGD ^[26] is typically indicator-based algorithm for solving MaOPs. MaOEA-IT ^[46] is a recently proposed MaOEAs which employs two-stage strategy to take the convergence and diversity of solutions into account.

3.2.   Performance metrics

We select two widely used performance metrics to evaluate the performance of DAEA and other compared algorithms on MaOPs. More details are as follows:

(1) Inverted generational distance (IGD) ^[57]: The IGD metric can measure both the convergence and diversity of the obtained solution set, whose mathematical description is:

where $P^*$ is a set of solutions obtained by uniformly sampling from the true PF and $A$ refers to non-dominated solutions obtained by an algorithm. $|P^*|$ is the size of $P^*$, and $d({f_i}^*, f)$ represents the Euclidean distance from the point in true PF to their closest to the obtained population $A$. A smaller IGD value of a set $A$ means a better performance of the set $A$.

(2) Hypervolume (HV) ^[27]: HV is a comprehensive indicator, which can assess both the convergence and diversity of a solution set as well. Its mathematical description is as follows:

where VOL$(\cdot)$ represents the Lebesgue measure, $A$ is the obtained non-dominated solutions, and $z = (z_1, \cdots, z_m)$ is the reference point and set to (1, 1, $\cdot$, 1). In addition, When the number of objectives $m$ exceeds 8, we employ the Monte Carlo method with 1, 000, 000 sampling points to approximate the HV values ^[25]. Obviously, the larger the HV value, the higher the quality of solution set $A$ is.

3.3.   Experimental settings

In this section, we outline the general experimental settings and provide specific parameter settings for each algorithm, as described below.

(1) Population size: The population size $N$ of these algorithms using reference vectors (such as NSGA-III, $\theta$-DEA) is determined by the parameter $H$ and the number of objectives $m$ ^[47]. To have a fair comparison, the population size of other compared algorithms keep identical to the above algorithms. Table 1 lists the detailed settings of population size $N$ for problems with different number of objectives used in this paper.

(2) Number of runs and termination criterion: Each algorithm is executed independently for 20 runs on every test instance, and the maximum function evaluations (MFEs) is used as the terminal criterion of each run. More details are presented in Table 1.

(3) Significance test: To evaluate the statistical significance of the results obtained from DAEA and the other five compared algorithms, the Wilcoxon rank sum test is utilized at a significance level of 0.05. In this test, we use the symbols "+", "-" and "=" to signify that the results obtained by a compared algorithm are significantly better than, worse than, and similar to 3DEA, respectively.

(4) Parameters of genetic operators: SBX ^[50] and PM ^[51] are used in all algorithms. For SBX ^[50], the crossover probability $p_c$ = 1.0 and the distribution index is set to $\eta_c$ = 20 are used in all compared algorithms. For the polynomial mutation ^[51], the distribution index and the mutation probability are set to $\eta_m$ = 20 and $p_m$ = $1$ for all compared algorithms.

(5) Parameter setting in each algorithm: The parameters used in all the compared algorithms remain unchanged from their original publications.

4.   Experimental results and analysis

4.1.   Performance comparisons on DTLZ test suite

Tables 2 and 3 present the statistical results of DAEA and five compared algorithms on the DTLZ test suite in terms of IGD and HV values, respectively, where the best results are highlighted in bold font.

From the IGD results in Table 2, it can be observed that DAEA is significantly superior to NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT on 26, 18, 35, 26 and 35 test instances, respectively. Analogously, it can be observed from the HV results in Table 3 that DAEA defeats NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT on 28, 18, 32, 20 and 35 test instances, respectively. In summary, our proposed 3DEA is superior to the other state-of-the-art algorithms.

Specifically, DTLZ1 is a multi-modal problem with the linear PF, which mainly challenges the algorithm ability to jump out of local optima. DAEA obtains the best results in dealing with this problem, followed by $\theta$-DEA.

For DTLZ2 test problem, DAEA performs best among the six compared algorithms. To more intuitively show the results obtained by six algorithms, we plot the parallel coordinates of final solutions obtained by each algorithm on 20-objective DTLZ2 instance in Figure 4. From this figure, we can see that $\theta$-DEA and DAEA have better comprehensive performance in terms of convergence and diversity, while the other algorithms fail to converge to the true PF.

Since its multimodal property, DTLZ3 is difficult to solve. From the IGD and HV values highlighted in Tables 2 and 3, we can see that our DAEA performs best for most instances with different number of objective except 16-objective instance, which are won by $\theta$-DEA. Other four algorithms are defeated by DAEA. It may because the double association in DAEA maintains a good diversity of solutions.

DTLZ4 mainly test the ability of algorithms to maintain population diversity. From Tables 2 and 3 we can see that DAEA performs best for most instances with different number of objective except 8-objective instance in Tables 3 and 8-, 12-, 16-objective instances in Tables 2, which are won by MaOEA-IGD. Other four algorithms are defeated by DAEA.

DTLZ5 and DTLZ6 have the degenerated PFs, which are used to test the ability of algorithms to deal with MaOPs with degenerate PFs. On these two problems, the IGD values obtained by DAEA are best for all instances. For HV metric, MaOEA-IGD is more suitable for solving this kind of problems than other compared algorithms. DAEA is the second and defeats the other four competitors on all test instances, which may be caused by the effective fitness assignment proposed by DAEA, since their tacit cooperation enables the algorithm to balance convergence and diversity.

DTLZ7, as a kind of mixed (concave and convex), disconnected and multi-modal problems, is difficult to better balance the convergence and diversity for existing evolutionary algorithms. From Tables 2, we can see that NSGA-III achieves the best IGD result on 5-objective instance, $\theta$-DEA achieves the best IGD results on 8-, and 12-objective instances, MaOEA-IGD has the best IGD statistical results on 16-, and 20-objective instances. DAEA is statistically similar to that obtained by $\theta$-DEA. From the HV results in Table 3, it can be observed that NSGA-III performs best on 5-objective instance, while $\theta$-DEA performing best HV on 8-, and 20-objective instances. For our DAEA, it performs best on 12-, and 16-objective instances. As for MaOEA-IT, it is defeated by its peer competitors.

4.2.   Performance comparisons on WFG test suite

As evidenced by statistical results of the IGD and HV values summarized in Tables 4 and 5, it is clearly observed that our DAEA performs best compared with other five algorithms in most cases. According to the statistic results of Wilcoxon rank sum test in the last line of Table 4, we can see that, among the 45 test instances, the IGD results obtained by DAEA are significantly better than those obtained by NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT in most test instances (in 25, 19, 34, 39 and 45 test instances, respectively), while the IGD results obtained by NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT are significantly better than those obtained by DAEA in only a few test instances (in 10, 9, 7, 3 and 0 test instances, respectively). According to the statistic results in Table 5, it can be concluded that DAEA defeats its competitors on most test instances. To be specific, DAEA outperforms NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT on 25, 21, 22, 39 and 42 test instances out of 45 test instances, respectively. Moreover, NSGA-III, $\theta$-DEA, SPEA/R, MaOEA-IGD and MaOEA-IT outperform DAEA in terms of HV results on a limited number of test instances (8, 6, 17, 1 and 0 test instances, respectively). In what follows, we will delve into detailed discussions regarding the experimental results.

WFG1 has a flat bias and a combination of convex and concave structures in its PF. From Table 4, we can observe that DAEA obtains the best IGD values on 5-objective instance, $\theta$-DEA achieves the best performance on 8-, 12- and 16-objective instances in terms of IGD values, while NSGA-III wins on 20-objective instance. In terms of the HV metric, SPEA/R exhibits the best performance on 5- and 8-objective instances, $\theta$-DEA shows the best performance on 12- and 16-objective instances, while NSGA-III demonstrates the best performance on the 20-objective instance. As for MaOEA-IGD and MaOEA-IT, both of them perform worse for this problem.

WFG2 consists of multiple convex PF segments that are not connected, and its variables cannot be separated.For this problem, NSGA-III and SPEA/R are superior to other four algorithms, followed by DAEA and $\theta$-DEA.

WFG3 is a variant of WFG2 where its PF is both linear and degenerate, which poses a big challenge to reference vector-based algorithms. Our proposed DAEA outperforms other methods on 8-, 12-, 16- and 20-objective test instances, while the best performance on 5-objective instance is achieved by $\theta$-DEA. The results verifies the effectiveness of DAEA on degenerate problems.

The problems WFG4 to WFG9 share the same hyper-ellipse PF shape in the objective space, but their characteristics differ in the decision space. To be specific, WFG4 has the property of multimodality with "hill sizes", it can be used to test whether a algorithm has the ability to jump out of local optima.

DAEA shows distinct advantage in most instances except the 8-objective instance, which is defeated by SPEA/R. WFG5 is a deceptive problem, the IGD metric shows that DAEA performs best on most test instances, while NSGA-III and $\theta$-DEA perform similar to DAEA on 8-, 12- and 16-objective instances. The HV metric shows that SPEA/R and DAEA are superior to the other four algorithms. WFG6 is a nonseparable and reduced problem, DAEA obtains the best IGD values on all test instances except 12-objective instance which is obtained by NSGA-III, which indicates that our proposed DAEA has promising performance of diversity and convergence, especially for the optimization problem with a larger number of objectives. WFG7-WFG9 all introduce some biases to challenge algorithms' diversity, DAEA performs best on these problems, followed by SPEA/R, $\theta$-DEA and NSGA-III, MaOEA-IGD and MaOEA-IT have poor performance for these problems.

4.3.   The effectiveness of double association

To verify the advantage of our double association, we use it to replace the association method of other decomposition-based algorithms. Here, $\theta$-DEA ^[16] and RPD-NSGA-II ^[17] as two classical decomposition-based MaOEAs are selected to make experiments. Both of them use single association method to associate solutions.

Tables 6 and 7 show the statistical results of $\theta$-DEA with different association methods in terms of IGD and HV values on DTLZ test suite, respectively. $\theta$-DEA_rand and $\theta$-DEA_DA represent using the random association method and double association method to replace the association method used in $\theta$-DEA. The best results are highlighted in bold. We can clearly observe from Tables 6 and 7 that after using our proposed double association method to replace the single association method used in $\theta$-DEA, the performance of $\theta$-DEA is significantly improved, while the random association makes the performance of $\theta$-DEA worse.

Tables 8 and 9 show the statistical results of RPD-NSGA-II with different association methods in terms of IGD and HV values on DTLZ test suite, respectively. RPD-NSGA-II_rand and RPD-NSGA-II_DA represent that using random association method and double association method to replace the association method used in RPD-NSGA-II. The best results are highlighted in bold. It can be seen clearly from Tables 8 and 9 that after using our proposed double association method to replace the single association method used in RPD-NSGA-II, the performance of RPD-NSGA-II is significantly improved, while the random association makes the performance of RPD-NSGA-II worse.

These results indicate that double association, as a new method to associate the solutions in objective space, can better maintain the diversity of the obtained solutions and give some potential solutions an opportunity for further evolution, and thus can improve the performance of reference point-based algorithms.

4.4.   Effect of different penalty parameter $\mu$

In our proposed quality evaluation scheme, a penalty parameter $\mu$ is designed to punish the solutions in dense area from global perspective. Considering that in an ideal state, each subspace should have two solutions distributed, so we set $\mu = \frac{|S(i)|}{2}$ to balance the convergence and diversity. To investigate the effect of different value of parameter $\mu$ on the proposed quality evaluation scheme, DAEA is tested on DTLZ1-7 with three different $\mu$ values: 1)$\mu = 1$; 2)$\mu = 5$; 3)$\mu = \frac{|S(i)|}{2}$.

Tables 10 and 11 show the statistical results of DAEA with three different $\mu$ values in terms of IGD and HV values on DTLZ test suite, respectively. The best results are highlighted in bold. As can be observed, DAEA_$\mu = \frac{|S(i)|}{2}$ has achieved most of the best results out of 35 test instances. From the statistic results of Wilcoxon rank sum test in the last line of Tables 10 and 11, we can also see that, among the 35 test instances, the IGD results obtained by DAEA_$\mu = \frac{|S(i)|}{2}$ are significantly better than those obtained by DAEA_$\mu = 1$ and DAEA_$\mu = 5$ in most test instances. These experimental results indicate that setting the penalty parameter $\mu$ as a dynamic parameter is much better than taking constants 1 and 5.

4.5.   Further discussion on DAEA

According to the experimental results (the effectiveness of double association) and comparative results, we can clearly see that these algorithms adopting the existing association operator for many-objective optimization may not work effectively and efficiently, especially in solving these problems with the multi-modal property or the complicated PF. One possible reason is that these algorithms ignore the importance of the subspace without associated, and further deteriorate their ability to explore unknown regions. In DAEA, we design the double association strategy that considers the empty subspace and associates it with a promising solution. This is helpful for DAEA to explore more unknown regions, and further improve the performance of DAEA. The effectiveness of the double association strategy has been validated by the corresponding experiment (i.e., the effectiveness of double association). In addition, most existing algorithms only consider the convergence and global diversity to measure the quality of solutions, which is not conductive to trade off the convergence and diversity of population. On the contrary, DAEA takes the convergence, global diversity, and local diversity into account, while using the designed dynamic penalty coefficient to adjust these three factors, and further balance the convergence and diversity of population. The effectiveness of the dynamic penalty coefficient has been confirmed by the corresponding experiment (i.e., parameter analysis).

However, our experimental results on the DTLZ7 indicate that the proposed DAEA does not present a clear advantage over $\theta$-DEA and MaOEA-IGD. Even so, the proposed DAEA still outperforms the other competitors. In addition, when DAEA deals with WFG1 and WFG2 with irregular PFs, DAEA is slightly worse than NSGA-III and SPEA/R. This is because of this fact that, for many-objective problems with the irregular Pareto fronts, the even distribution of reference vectors combing with our proposed double association strategy may result in partial waste of computing resources and affect the performance of algorithm.

5.   Conclusions and future work

In this paper, we have presented a double association-based evolutionary algorithm (denoted as DAEA) for many-objective optimization problems. The proposed double association strategy takes the empty subspace into account and associates it with a solution that is closest to this subspace, which can increase the probability of an unknown area being explored. In addition, a new quality evaluation scheme that takes the convergence and diversity of solutions into account is developed to measure the quality of each solution in subspace, where the diversity of each solution is subdivided into global diversity and local diversity. Then, a dynamic penalty coefficient is designed to protect these solutions located in sparse areas by penalizing the worse global diversity of solutions.

To demonstrate the high competitiveness of DAEA, we compare it with five state-of-the-art many-objective evolutionary algorithms on two test suites, DTLZ and WFG, with the number of objectives varying from 5 to 20. The corresponding experimental results demonstrate that our proposed DAEA has higher competitiveness in terms of both convergence enhancement and diversity maintenance compared with the other state-of-the-art MaOEAs.

Although DAEA outperforms the compared algorithms on most of the test problems, it struggles in dealing with some problems with the irregular Pareto fronts. For future research, we would like to add the idea of reference vector dynamic adjustment mechanism on DAEA to improve the algorithm's adjustability on the irregular Pareto fronts. In addition, we also consider verifying the performance of DAEA on MaOPs in the real world.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by National Natural Science Foundation of China (NO. 62202366) and Natural Science Basic Research Program of Shaanxi Province of China (2022JQ-624 and 2021JQ-674).

Conflicts of interest

The authors declare that they have no conflicts of interest.