A Frequency-based Parent Selection for Reducing the Effect of Evaluation Time Bias in Asynchronous Parallel Multi-objective Evolutionary Algorithms

Parallel evolutionary algorithms (PEAs) have been studied for reducing the execution time of evolutionary algorithms by utilizing parallel computing. An asynchronous PEA (APEA) is a scheme of PEAs that increases computational efficiency by generating a new solution immediately after a solution evaluation completes without the idling time of computing nodes. However, because APEA gives more search opportunities to solutions with shorter evaluation times, the evaluation time bias of solutions negatively affects the search performance. To overcome this drawback, this paper proposes a new parent selection method to reduce the effect of evaluation time bias in APEAs. The proposed method considers the search frequency of solutions and selects the parent solutions so that the search progress in the population is uniform regardless of the evaluation time bias. This paper conducts experiments on multi-objective optimization problems that simulate the evaluation time bias. The experiments use NSGA-III, a well-known multi-objective evolutionary algorithm, and compare the proposed method with the conventional synchronous/asynchronous parallelization. The experimental results reveal that the proposed method can reduce the effect of the evaluation time bias while reducing the computing time of the parallel NSGA-III.


Introduction
Evolutionary algorithms (EAs) have been applied to a wide range of real-world optimization problems owing to their high search capability without any problem-specific knowledge. When applying EAs to real-world applications, solution evaluations may take much computational time, such as due to physical simulation or complex numerical calculations.
Parallel EAs (PEAs) [1][2][3] are a promising technique to speed up the optimization process for computationally expensive problems. A masterworker parallelization [4] is one of the typical approaches of PEAs, where a single master computing node executes the main procedure of EA, e.g., initialization, parent selection, genetic operations, and survival selection. In contrast, many worker nodes evaluate each solution in parallel.
Master-worker PEAs (MW-PEAs) can be mainly classified into two approaches, a synchronous PEA (SPEA) and an asynchronous PEA (APEA) [5]. SPEA generates a population in the next generation after evaluating all solutions. On the other hand, in APEA, a new solution to be evaluated is generated immediately after completing one solution evaluation. Since SPEA needs to wait for a solution with the longest evaluation time before generating the next population, it increases the idling time of worker nodes and decreases the computational efficiency. On the other hand, APEA can overcome this drawback because it can continuously evolve solutions without the idling time of worker nodes. However, since APEA generates a new solution whenever a solution evaluation completes, it could lead to local optima with a short evaluation time [6,7]. This paper proposes a new parent selection method to reduce the effect of evaluation time bias in APEA. Concretely, the proposed method considers the search frequency of each search region and selects parents so that the search frequency of all solutions becomes uniform. The proposed method introduces a new parameter that stores the search frequency of solutions and selects parents from solutions with fewer search frequencies.
This paper is an extended version of the author's work [8]. This paper improves the behavior analysis of the proposed method in more detail, mainly: the previous work analyzed the behavior of the proposed method using only one multi-objective benchmark problem (DTLZ1). On the other hand, to achieve a deeper analysis, this study utilizes multimodal multi-objective test functions (MMFs) [9] and designs benchmarks with two Pareto-optimal solution sets where each Pareto-optimal solution has a different evaluation time. Comparing the proposed method with synchronous/asynchronous parallel MOEA shows that the proposed method can equally obtain Pareto-optimal solution sets with different evaluation times (see Section 4). In addition, this paper provides further analysis of a proposed method parameter, mainly r s in the proposed method (see Section 3.2).
The rest of this paper is organized as follows. The following section briefly introduces APEAs and mentions their problems in the evaluation time bias. Section 3 proposes the parent selection strategy and shows its concrete example on NSGA-III. Section 4 defines the test problems used in this work, and Section 5 describes the experimental settings. Then, the parameter setting of the proposed method is discussed in Section 6, and Section 7 compares the proposed method with SPEA and APEA. Finally, Section 8 concludes this paper and addresses future works.

Background
Parallel evolutionary algorithms (PEAs) [1][2][3] have been studied to reduce the computing time of EA methods by executing a single EA run on multiple computing nodes and have been applied to several real-world applications, such as education [10], data mining [11], nanoscience [12], and routing [13].
A master-worker model (known as a global model) is a straightforward approach to implementing PEAs [4,14], and is widely used in many recent works [15,16]. On a master-worker PEA (MW-PEA), a master computing node executes the main procedure of EAs, such as initialization, selection, genetic operators, and replacement. On the other hand, many worker nodes evaluate newly generated solutions in parallel and return their results to the master node.
MW-PEAs can be classified into synchronous PEAs (SPEAs) and asynchronous PEAs (APEAs). SPEA waits for all evaluations of solutions executed by worker nodes and generates a new population using all newly evaluated solutions. Since SPEA needs to wait for the longest evaluation for each generation, computational efficiency decreases if the evaluation times differ. On the other hand, APEA continuously generates a new solution without waiting for evaluations of other solutions. This enables the efficient use of the computing resource even if the evaluation times of solutions differ. Since previous research has demonstrated the effectiveness of APEAs, for example, on continuous optimizations [17], genetic programmings [18], and multi-objective optimizations [19], this paper focuses on APEAs.
Although APEA is a practical approach of MW-PEAs, the previous studies demonstrated that APEA is biased toward the search region having a short evaluation time if the evaluation time differs depending on the search region [6,7]. This happens because APEA gives many opportunities to search for solutions with a short evaluation time. SPEA, on the other hand, is an option to avoid the effect of the evaluation time bias because it is not affected by the evaluation time bias. However, SPEA still wastes waiting time because the evaluation times of solutions differ.
From the above, it can be seen that SPEA has poor computational efficiency regardless of the evaluation time bias. On the other hand, APEA is computationally efficient but is affected by the evaluation time bias. Therefore, this research proposes a method to reduce the effect of the evaluation time bias in APEA while maintaining its computational efficiency.

Proposed method
This section proposes a new parent selection method for reducing the effect of evaluation time bias in APEAs. The following subsection first explains the basic concept of the proposed method, and then Section 3.2 introduces the proposed parent selection. Finally, Section 3.3 shows an example of applying the proposed method to the asynchronous parallel NSGA-III.

Basic concept
APEAs are negatively affected by the evaluation time bias because the search frequency for solutions with a short evaluation time increases. On the other hand, SPEAs are not affected by the evaluation time bias because the search frequency is almost the same in all search regions regardless of the evaluation time bias. This fact suggests that adjusting the search progress of all solutions to be uniform can effectively reduce the effect of evaluation time bias in APEAs.
From this viewpoint, this paper proposes a new parent selection method that introduces a new parameter, a search frequency parameter, that stores how many offspring are generated from each solution. The proposed selection method selects parents to preserve the uniformity of the search frequency of solutions according to the additional frequency parameter. This contributes to preventing excessive parent selection of solutions in the regions with a short evaluation time. Fig. 1 illustrates the conventional and proposed APEAs on a one-dimensional maximization problem with the evaluation time bias. This example has a global optimum with a long evaluation   Fig. 1a increases the search frequency of solutions with a short evaluation time (the right area). This eliminates solutions close to the global optimum (a long evaluation time) by comparing them with the more frequently searched solutions near the local optimum. To overcome this problem, the proposed method depicted in Fig. 1b stores the search frequency of solutions (s.f req) and attempts to preserve the uniformity of these frequencies. This may avoid eliminating solutions with a long evaluation time due to its slow search progress.

Algorithm
The proposed method introduces a new parameter to store the search frequency of each solution. Let the search frequency of a solution s be s.f req. First, the frequency parameter s.f req for all solutions is initialized to 1. After that, when a solution is selected as a parent, the frequency parameter is incremented as s.f req ← s.f req + 1. In addition, generated offspring s new inherits the frequency parameter as the mean of its parents. In particular, when an offspring solution s new is generated from two parents, p 1 and p 2 , the search frequency parameter of s new is calculated as s new .f req ← (p 1 .f req + p 2 .f req)/2. This allows us to store the Algorithm 1 A pseudo-code of a simple APEA with the proposed method. The underlined texts are specific to the proposed method.
1: Generate S random solutions ⊲ S is the number of worker nodes 2: Send solutions to worker nodes 3: P 0 ← ∅ 4: while |P 0 | < N do ⊲ N is the population size 5: s ←wait for a solution from worker nodes 6: Generate a random solution s 8: s.f req ← 1

9:
Send s to an idling worker node 10: end while 11: t ← 0 12: while Terminal conditions do 13: s c ←wait for a solution from worker nodes 14: P t+1 ← select N solutions from P t ∪ {s c } ⊲ Any replacement strategy 15: Sort P t+1 in ascending order according to s.f req (s ∈ P t+1 ) 16:

21:
Send s new to the idling worker node 22: t ← t + 1 23: end while search progress of each search region as additional information for each solution.
Unlike the conventional parent selection, the proposed method selects only solutions with fewer search frequencies. In particular, when the parent selection, the proposed method preliminary selects a candidate pool from the current population according to the search frequency parameter. This can reduce the search opportunities for solutions with a short evaluation time and lead solutions to be uniformly selected as the parents.
Algorithm 1 shows a pseudo-code of an APEA with the proposed parent selection. Additional procedures from a simple APEA are underlined. The search frequency parameters of the initial solutions are set to 1 (Step 8). The population is sorted in ascending order of the frequency parameters when selecting parent solutions (Step 15). Then, the top r s |P t+1 | solutions in P t+1 are extracted as a parent candidate pool P ′ t+1 in Step 16. Here, r s (0 ≤ r s ≤ 1) is a selection ratio parameter that determines how the uniformity of search frequency is prioritized. Since this may affect the search capability of the proposed method, its effect will be discussed in Section 6. Then, parents are selected from P ′ t+1 according to the algorithm-specific selection method (e.g., tournament selection, roulette-wheel selection). After generating an offspring, the frequency parameters of the parents are incremented by one (Step 19), and a newly generated solution inherits the frequencies of the parents ( Step 20).
The difference between the conventional APEA and the proposed method is that: the conventional one selects parents from the entire population regardless of the search frequency. This induces that solutions having short evaluation times frequently get opportunities to be selected as parents, and the search direction is biased. On the other hand, the proposed method considers the search frequency of solutions and selects parents from less selected solutions for the offspring generation. This mechanism can allow selection as parents for all solutions and prevent the asynchronous evolution from being affected by the evaluation time bias.

An example of the proposed method with NSGA-III
NSGA-III [20] is one of the most well-known and successful MOEA methods combining dominance and decomposition strategies. This section shows an example of applying the proposed method to the asynchronous parallel NSGA-III, which will be used in experiments in Sections 6 and 7. Please see the detailed algorithm in [20]. Algorithm 2 describes the brief flow of the asynchronous parallel NSGA-III assisted by the frequency-based selection (FS-NSGA-III). The master node initializes the population (Step 2) and assigns the search frequency parameter s.f req = 1 for all solutions in the initial population P 0 (Step 3). After the initialization, the master node sends all solutions to worker nodes (Step 4), and the main procedure is repeated until satisfying the termination condition. When generating offspring, the proposed method sorts the population in ascending order of the search frequency parameter (Step 7) and selects the top Algorithm 2 An algorithm of the asynchronous parallel NSGA-III with the proposed frequencybased parent selection 3: s.f req = 1 (s ∈ P 0 ) 4: Send all solutions to worker nodes. 5: Wait for evaluations of all solutions. 6: while Termination conditions do 7: 10: : Send s new to an idling worker node. 14: s ← wait for the next evaluation.

15:
R t ← P t ∪ {s}. 16: 17: In contrast with NSGA-III randomly selecting two parent solutions from the entire population, FS-NSGA-III selects parents from the candidate pool P ′′ t (Step 9). After generating offspring, the proposed method increments the frequency parameters of the parents as p i .f req ← p i .f req + 1 (Step 11) and inherits the frequency parameter of the generated offspring as the mean of its parents as s new .f req ← (p 1 .f req + p 2 .f req)/2 (Step 12). Then, the master node sends the generated offspring to an idling worker node (Step 13) and waits for the subsequent evaluation of a solution (Step 14). When receiving an evaluation, NSGA-III selects the next population from the current population and a newly evaluated solution (Step 16). NSGA-III uses the front ranking and the niche-preservation operation based on the reference point in the selection procedure.

Test problems with evaluation time bias
This work designs multi-objective optimization test problems with the evaluation time bias to deeply analyze the behavior of the proposed method. In particular, this work uses multimodal   multi-objective test functions (MMFs) [9], which are bi-objective optimization problems with multiple separate Pareto sets (PS) in different regions of the decision space. One of the notable features of MMFs is that the Pareto front in the objective function space can be entirely approximated by only acquiring one of the separated Pareto sets. This study uses two-dimensional MMF2-6 and MMF8 1 , which have two separated Pareto sets. Fig. 2 depicts the Pareto set for each problem (see [9] for more detailed problem definitions). The Pareto set with a smaller x 2 value for each x 1 is defined as PS1, while another one is defined as PS2 as follows: where x 1 and x 2 are the first and second design variables of x, while θ(x 1 ) is a problem-dependent function that returns a boundary plane calculated by x 1 . The functions c s (x) and c l (x) are defined Table 1: The definition of functions c s and c l , and a parameter σ for each problem as shown in Table 1, which is determined from the definition of the Pareto set.
Based on the Pareto set in MMFs, this work designs the evaluation time function so that the optimal solutions in PS2 require a longer evaluation time than those in PS1. Specifically, the biased evaluation time is defined as: This work names t bias Bias. For each problem, σ determines the variance of the Gaussian function, and the value of σ is shown in Table 1. An example of the evaluation time distribution of MMF2 is shown in Fig. 3. The horizontal axis represents x 1 , while the vertical axis represents x 2 . The color bar indicates the evaluation time. In this setting, solutions in the PS1 have a shorter evaluation time than those in the PS2. It is expected that APEA will converge more quickly to the PS1 in such a situation.
In addition to Bias, this experiment uses a non-biased evaluation time function that returns a random value sampled from the normal distribution as t norm (x) ∼ N (t p , c v t p ). The variable t p denotes the mean evaluation time, while c v determines the variance of the evaluation time. Since t norm (x) is independent of the decision variable,

Experimental setting
This paper conducts experiments on the simulated parallel computational environment to investigate the effectiveness of the proposed method. The proposed method is applied to NSGA-III as a concrete algorithm shown in Section 3.3. The experiments compare three parallel NSGA-IIIs, synchronous parallelization (SP-NSGA-III), asynchronous parallelization (AP-NSGA-III), and the proposed parallelization (FS-NSGA-III). Note that this paper does not aim to solve multimodal multi-objective optimization problems efficiently but uses them to analyze the influence of evaluation time bias. Thus, this study does not uses specific techniques for finding niches.
This section first explains the simulation environment used in the experiment, and then the parameter settings used in the experiment are shown. The final subsection provides evaluation criteria for assessing the competitive methods.

Simulated parallel computational environment
The experiments use a simulated parallel computational environment based on the computational time model proposed in the work of [21]. This model consists of a single master node and λ worker nodes. The master node computes the main procedure of the EA algorithm in t s = 1 simulation time. In contrast, the worker nodes evaluate one solution and return their evaluation results. This experiment simulates λ = 100 worker nodes where 100 solutions are simultaneously evaluated. The evaluation times on the worker nodes depend on the Bias and No-bias functions. In Bias, the value of t mean is set to 1000. In such a setting, the maximum evaluation time is almost 2000, while the minimum one is almost one, so the longest evaluation time is 2000 times longer than the shortest one. On the other hand, in No-bias, the mean evaluation time t p = 1000, while the variance parameter c v = 0.2.

Parameters
The experiments were conducted for 31 independent runs for each parallelization method.

Evaluation criteria
This experiment uses the inverted generational distance (IGD) indicator [22] to assess the quality of the obtained solutions in the objective space.
The IGD value is calculated as: where P * denotes a reference point set (the true Pareto solution set), while P denotes the nondominated solutions obtained by the algorithm. The function d(x, y) calculates the Euclidean distance between x and y. The solutions obtained by the algorithm are worthful if the IGD value is small.
In addition, the IGDX indicator [23] is used to evaluate the quality of solutions in the design variable space. The IGDX value is calculated as: where P * and P denote the true Pareto solution set and the obtained one by the algorithm. When calculating IGD, the distance is calculated on the objective space. On the other hand, when calculating IGDX, the distance on the design variable is calculated. This work independently calculates the IGDX values for two separate Pareto sets to confirm whether both Pareto sets are obtained simultaneously. The IGDX values for PS1 and PS2 are denoted as IGDX 1 and IGDX 2 , respectively. To assess if both Pareto sets are equally obtained, the difference between two IGDX values is defined as: ∆IGDX(P * , P ) = IGDX(P S1, P ) − IGDX(P S2, P ). (4) If ∆IGDX = 0, both Pareto sets are equally obtained. On the other hand, if ∆IGDX < 0, since IGDX 1 < IGDX 2 , the algorithm is biased to PS1, and vice versa. In the experiment using Bias, the ∆IGDX value of SP-NSGA-III is expected to be 0, while that of AP-NSGA-III may be less than 0 because its search direction is biased to PS1. It can be expected that the proposed method shows similar behavior to the synchronous method by reducing the effect of the evaluation time bias.
The Kruskal-Wallis test will be performed to confirm a statistical difference between the three parallelization methods for each criterion. If a significant difference is found with the Kruskal-Wallis test, we perform the post-hoc pairwise comparisons using the Wilcoxon rank-sum test with the Bonferroni adjustments. The horizontal axis shows the selection ratio r s , while the vertical axis shows the IGD value. The two boxes are connected with the "*" symbol if a significant difference with a significance level of 5% is found, and the "**" symbol if the significance level is 1%. These results show that r s = 0.1 obtains a significantly worse IGD value in MMF4, MMF5, and MMF6. Moreover, the selection ratio of r s = 0.1 shows relatively worse IGD values than the other settings in other problems. These results indicate that a small r s negatively affects the search capability in No-bias and Bias. When using r s = 0.1, the proposed method selects parents from 10% of solutions in the current population, which restricts the diversity of the parents and reduces the search capability.
On the other hand, the results with r s ≥ 0.3 show no significant difference in all problems and both evaluation times. These results indicate that a small selection ratio should be avoided, but larger selection ratios do not essentially affect the search capability.  table summarizes the median and IQR values. The horizontal axis shows the selection ratio r s , while the vertical axis shows the ∆IGDX value. Like figures 4 and 5, the two boxes with a significant difference are connected with the "*" or "**" symbols.

Effect of evaluation time bias
First, the results with No-bias show that there is no significant difference between different selection ratios. On the other hand, Fig. 7 shows that the ∆IGDX value decreases as the selection ratio r s increases on Bias. Specifically, the selection ratio of r s = 0.9 is significantly biased toward the search region with a shorter evaluation time (a negative ∆IGDX value). Since a large selection ratio selects parent individuals from a large candidate pool containing solutions with a large search frequency, it gets close to the standard asynchronous method. Therefore, the larger the selection ratio r s , the more susceptible the evaluation time bias. These results indicate that a large selection ratio should be avoided to reduce the effect of the evaluation time bias.

Computational efficiency
This subsection analyzes the computational efficiency of different selection ratios by comparing the simulation time until obtaining a particular quality of the Pareto front. Concretely, the target IGD value is decided for each problem from the previous results as follows: 2.0 × 10 −3 in MMF2 and MMF3, 2.0 × 10 −4 in MMF4, 4.5 × 10 −4 in MMF5, 3.5 × 10 −4 in MMF6, and 1.5 × 10 −4 in MMF8. Figures 8 and 9 show the simulation time until the target IGD value is reached for No-bias and Bias, and the bottom tables summarize the median and IQR values of simulation time until reaching the target IGD value. The horizontal axis shows the selection ratio, while the vertical axis shows the simulation time. As with the previous results, two boxes with a significant difference are connected with the "*" symbols.
First, no significant difference is found in all test problems when using No-bias. However, the selection ratio r s = 0.1 requires a relatively longer simulation time than the others in all problems. For the other selection ratios of r s ≥ 0.3, all selection ratios take almost equal simulation time in No-bias. These results can be explained because a small selection ratio of r s = 0.1 shows the lower search capability, as demonstrated in Section 6.1.
From the results with Bias, on the other hand, the larger the selection ratio is used, the shorter the simulation time is obtained. Here, it is necessary to consider the effect of the evaluation time bias on the execution time. Specifically, when using the selection ratios of r s = 0.9, the search direction is biased toward a region with a shorter evaluation time, as demonstrated in Section 6.2, resulting in a shorter execution time. In fact, it can be seen that the evaluation time of the selection ratios of r s = 0.9 is the shortest in all problems, and some are significantly shorter than others. For the other selection ratios, the selection ratio of r s = 0.1 obtains the significantly longest execution time in MMF4, MMF5, and MMF6. This is due to a combination of the fact that the small selection ratio decreases its search capability (shown in Section 6.1), and it obtains solutions with both longer and shorter evaluation times by reducing the effect of the evaluation time bias (shown in Section 6.2). In contrast, the selection ratios of r s = 0.5, 0.7 acquire stably shorter execution times in both evaluation times. These ratios acquire comparable search capability to the large selection ratios (e.g., r s = 0.9) and decrease the influence of the evaluation time bias by using a smaller ratio (e.g., r s = 0.1). For this fact, it is indicated that such moderate selection ratios are appropriate.

Comparison of Different Parallelization Schemes
This section compares the performance of different parallelization schemes, SP-NSGA-III, AP-NSGA-III, and FS-NSGA-III. Since the results in the previous section suggested the selection ratio such as r s = 0.5, 0.7, the following experiments use r s = 0.5 in the proposed method. Figures 10 and 11 show the boxplot of the IGD value after the maximum number of evaluations for No-bias and Bias, and the bottom tables summarize the median and IQR values. The horizontal axis shows the different parallelization methods, while the vertical axis shows the IGD value. As in the previous section, a significant difference is depicted with "*" symbols. From these results, when using No-bias, there is no significant difference between the three parallelization methods. On the other hand, from Fig. 11, no significant difference is found except for MMF8 when using Bias. FS-NSGA-III produces a comparable IGD value in other problems compared with other methods. In MMF8, FS-NSGA-III obtains a significantly larger (worse) IGD value than AP-NSGA-III. This result can be explained in Fig. 13. Since AP-NSGA-III is biased toward the search region with a short evaluation time, it precisely obtains the Pareto front by only solutions in PS1. In contrast, since FS-NSGA-III and SP-NSGA-III obtain both Pareto sets equally, the IGD values are inferior to those of AP-NSGA-III that approximates the Pareto front elaborated by PS1 only.  This result indicates that FS-NSGA-III does not negatively affect the search capability of AP-NSGA-III, even though selecting the parents from the limited candidate pool. Figures 12 and 13 show the boxplot of the ∆IGDX value after the maximum number of fitness evaluations for No-bias and Bias, and the bottom tables summarize the median and IQR values. The horizontal axis shows the different methods, while the vertical axis shows the ∆IGDX value. As in the previous section, a significant difference is depicted with "*" symbols.

Effect of evaluation time bias
First, the No-bias results show no significant difference in the ∆IGDX value between the three parallelization methods. In addition, since the ∆IGDX value is almost zero in all problems, it is revealed that all parallelization schemes can obtain the separated Pareto sets equally if the evaluation time is not biased.
On the other hand, when using Bias, significant differences are found in MMF2, MMF4, MMF5, and MMF8, while no significant difference is found in MMF3 and MMF6. In particular, AP-NSGA-III obtains a significantly smaller (negative) ∆IGDX value than SP-NSGA-III in MMF2, MMF4, MMF5, and MMF8. This brings out the effect of the evaluation time bias in the asynchronous method.
From the results of the proposed method, a significant difference in MMF4 and MMF5 can be found. In these problems, the ∆IGDX value of FS-NSGA-III is not significantly different from that of SP-NSGA-III. In contrast, AP-NSGA-III is significantly biased toward PS1 (shorter evaluation time) than FS-NSGA-III and SP-NSGA-III. On the other hand, in MMF2 and MMF8, FS-NSGA-III shows a significantly smaller (negative) ∆IGDX value than SP-NSGA-III. Moreover, no difference between FS-NSGA-III and AP-NSGA-III is found, though the distribution of the ∆IGDX value of FS-NSGA-III is close to zero compared with AP-NSGA-III.
These results can be classified into three categories that are; (1) MMF3 and MMF6, where the ∆IGDX values are almost equal between the three methods (Proposed ≈ Sync. ≈ Async.); (2) MMF4 and MMF5, where the asynchronous method is significantly biased toward the region with a shorter evaluation time than the others (Proposed ≈ Sync. ≫ Async.); and (3)   and MMF8, where the proposed method is also biased (Sync. ≫ Proposed ≈ Async.). The difference between these categories can be explained from the perspective of the distribution of the Pareto set shown in Fig. 2. In MMF3 and MMF6, two Pareto sets are overlapped in the x 2 dimension, and a solution in one Pareto set is easily generated from a solution in another Pareto set. Thus, the evaluation time bias effect is small, and all methods are not biased. On the other hand, two Pareto sets are completely separated in the x 2 dimension in the other benchmarks, but they are close in MMF4 and MMF5 compared with MMF2 and MMF8. When the regions of Pareto sets are separated and their evaluation times are biased, the asynchronous method results in a biased search toward regions with short evaluation times. On the other hand, the proposed method can reduce the effect of evaluation time bias even when optimal solutions exist in separate regions with the biased evaluation time.
These results indicate that the proposed method can reduce the effect of the evaluation time bias despite being asynchronous. In contrast, the asynchronous method without the proposed method easily converges to a Pareto set with a shorter evaluation time.

Computational efficiency
This subsection analyzes the computational efficiency of three methods by comparing the simulation time until obtaining the target IGD values defined in Section 6.3. Figures 14 and 15 show the simulation time until obtaining the target IGD value for No-bias and Bias, and the median and IQR values are summarized at the bottom of the figures. The horizontal axis shows the different parallelization methods, while the vertical axis shows the simulation time. As in the previous section, a significant difference is depicted with "*" symbols.
From these results, FS-NSGA-III significantly reduces the execution time compared with the SP-NSGA-III when using Bias. Meanwhile, there is no significant difference between FS-NSGA-III and AP-NSGA-III. This indicates that the proposed method retains the computational efficiency of the asynchronous one.
On the other hand, from the results with the evaluation time of Bias, FS-NSGA-III also  acquires a shorter execution time than SP-NSGA-III. In particular, the proposed method significantly reduces the execution time in MMF4, MMF5, and MMF6. Although no significant difference is found in MMF2, MMF3, and MMF8, the proposed method obtains enough better performance in half execution time than SP-NSGA-III. Since the result in Fig. 13 showed that FS-NSGA-III and SP-NSGA-III equally obtain two Pareto sets with different evaluation times, it can be said that the proposed method can reduce the execution time while reducing the effect of evaluation time bias. The comparison of FS-NSGA-III and AP-NSGA-III indicates that the proposed method requires a significantly longer execution time when using Bias. However, this behavior can be explained because AP-NSGA-III is biased toward searching for solutions with shorter evaluation times (PS1), as indicated in Fig. 13. On the other hand, since FS-NSGA-III obtains solutions with longer evaluation times, its execution time increases compared with AP-NSGA-III. This result indicates that the proposed method obtains the computational efficiency of the asynchronous method while avoiding the effect of the evaluation time bias.

Conclusion
This paper proposed a new parent selection method for reducing the effect of evaluation time bias in APEAs. In particular, the proposed method considers the search frequency of each solution and selects parents from the preselected candidate pool. This paper conducted experiments on multi-objective optimization test problems based on MMFs to analyze the effect of the evaluation time bias deeply. The proposed method was applied to the parallel NSGA-III and was compared with the synchronous and the asynchronous parallelizations.
The experiments first analyzed the impact of the selection ratio in the proposed method using the same test problems. This analysis indicated that the selection ratio in 0.5 ≤ r s ≤ 0.7 acquires an appropriate balance between the search capability and the computational efficiency while reducing the effect of the evaluation time bias. Then, the experimental results indicated that the proposed method could reduce the negative influence of the evaluation time bias. The proposed method also does not adversely affect the search capability of APEAs while reducing the execution time significantly from SPEAs. These results revealed that the proposed method possesses high search capability and high computational efficiency for problems with heterogeneous evaluation time.
It should be addressed to further analyze the proposed method on other benchmarks and with other EA methods shortly. In addition, although this paper only compared the proposed method with the synchronous and the asynchronous method, it should be compared with or integrated into a semi-asynchronous method [24] to adapt any characteristics of the evaluation time.