A Reference Point-Based Evolutionary Algorithm Solves Multi and Many-Objective Optimization Problems: Method and Validation

The integration of a decision maker's preferences in evolutionary multi-objective optimization (EMO) has been a common research scope over the last decade. In the published literature, several preference-based evolutionary approaches have been proposed. The reference point-based non-dominated sorting genetic (R-NSGA-II) algorithm represents one of the well-known preference-based evolutionary approaches. This method mainly aims to find a set of the Pareto-optimal solutions in the region of interest (ROI) rather than obtaining the entire Pareto-optimal set. This approach uses Euclidean distance as a metric to calculate the distance between each candidate solution and the reference point. However, this metric may not produce desired solutions because the final minimal Euclidean distance value is unknown. Thus, determining whether the true Pareto-optimal solution is achieved at the end of optimization run becomes difficult. In this study, R-NSGA-II method is modified using the recently proposed simplified Karush–Kuhn–Tucker proximity measure (S-KKTPM) metric instead of the Euclidean distance metric, where S-KKTPM-based distance measure can predict the convergence behavior of a point from the Pareto-optimal front without prior knowledge of the optimum solution. Experimental results show that the algorithm proposed herein is highly competitive compared with several state-of-the-art preference-based EMO methods. Extensive experiments were conducted with 2 to 10 objectives on various standard problems. Results show the effectiveness of our algorithm in obtaining the preferred solutions in the ROI and its ability to control the size of each preferred region separately at the same time.


Introduction
Most real-world optimization problems usually contain two or more conficting objective functions. Tese objective functions must be optimized simultaneously. Tis type of problem is known as a multi-objective optimization problem (MOP).
In MOPs with contradictory objectives, a single solution that can be considered the best does not always exist. Instead, a set of solutions represents the best compromises among the diferent objectives. Tis set, which belongs to the search space, is known as the Pareto set (or efcient set), whereas its images, which belong to the objective space, are known as Pareto front (PF) [1,2]. Several evolutionary multi-objective optimization (EMO) algorithms, such as NSGA-II [3], SPEA2 [4], and MOEA/D [5], have been suggested in the past two decades or more. Classical EMO mainly aims to obtain a set of well-converged and well-distributed nondominated solutions that approach the entire PF [6,7]. Researchers have devoted their efort to developing algorithms in recent years [8][9][10][11][12][13][14].
Te proportion of non-dominant solutions rises as the number of objectives increases, which is one of the fundamental problems with all EMO approaches. Due to the insufcient selection pressure caused by a high percentage of non-dominant solutions, the EMO approach cannot advance in fnding the optimal spots. Incorporating the decision maker's preferences into the algorithm is a practical way to deal with this issue. A new kind of ranking mechanism [9] can be used to make selection pressure stronger and steer the optimization approach to search in a specifed region.
In real life, the DM is always focused on some specifed subsets of the obtained solutions. Te techniques of preference-based MCDM aim to fnd a part of the PF, whereas EMO algorithms aim to obtain a well-distributed set of points close to the whole PF. We call the part of the optimal solutions that is near to or lies on the PF a region of interest (ROI) [15]. Solutions within the ROI satisfy the DM's need. However, this scenario does not mean that any efcient solutions outside the ROI are not the optimal solutions to the problem. Te preference information given by the DM in the EMO can enable a highly efcient search. Many diferent approaches to preference information given by DMs exist, such as reference point (RP), preference angle, and reference weights. One of the most utilized approaches in preferencebased EMO algorithms is the RP. As mentioned above, EMO tries to fnd well-distributed multiple efcient solutions across the whole PF, as displayed in Figure 1(a). Also, this fgure illustrates the feasible objective region and the unfeasible objective region. On the contrary, preference-based EMO algorithms concentrate on a certain part of the true PF based on a preference point (reference point) determined by DM. Te non-dominated points cluster near the RP, as shown in Figure 1 Te following is a typical classifcation of methods based on preferences, depending on how they are expressed by the DM [16][17][18]: (i) a priori methods, where preferences are expressed before calculating PO solutions, for example, through a utility function [19] or by an RP [20]; (ii) a posteriori methods, in which the DM chooses the solution of her/his preference after a set of efcient solutions has been calculated (for example, [21,22]); (iii) interactive methods, where the DM guides the search with a utility function, and this function may change during the optimization process because of the new information acquired (for example, [23,24]); and fnally, (iv) methods not based on preferences, where additional information on preferences is not available, and the idea is to fnd a balance between the objectives [25].
Over the past two decades, researchers have focused their attention on preference-based EMO approaches. Tese approaches have been actively developed, and they mainly focus on specifc parts of the PF. Depending on the preference information supplied by the DM, these algorithms seek to fnd an ROI that is close to/on the true PF.
Numerous algorithms of preference-based EMO have been introduced. Deb and Sundar [26] suggested the RPbased NSGA-II (R-NSGA-II), which focuses on obtaining a preferred ROI during the evolutionary process. By including the RP's location information in the Pareto dominance, Molina et al. [27] initiated a concept of Pareto dominance termed g-dominance. Ben Said et al. [28] presented a novel variant of the Pareto-dominant relationship, called r-dominance, with which we can obtain good convergence to the PF. Ruiz et al. [29] proposed WASF-GA, another variant of the preference-based MOEA algorithm. Yu et al. [30] suggested a diferent representative preference-based decomposition MOEA by decomposing the preference information into several scalar optimization problems. Recently, new R-NSGA-II modifed methods have been proposed to assist DMs in convergent to Pareto-dominance compliant solutions in a specifc area of interest [31][32][33].
Although many preference-based algorithms use various metrics to select preferred solutions, some of these metrics require prior knowledge of the PF while others require specifc parameters [34,35]. S-KKTPM does not require prior knowledge of the PF or any parameters.
Herein, we introduce a novel preference-based NSGA-II algorithm. Te Euclidean distance was utilized in the original R-NSGA-II study as a metric between two trade-of solutions. In our study, we use the simplifed Karush-Kuhn-Tucker proximity metric (S-KKTPM) instead of the Euclidean distance metric. S-KKTPM can anticipate the convergence behavior of a point from PF without prior knowledge of the optimum solution [36,37]. Te Karush-Kuhn-Tucker (KKT) conditions occupy a signifcant role in optimization theory. KKT proximity measure was proposed through these conditions. Incorporating S-KKTPM within the R-NSGA-II provides theoretical convergence properties for the fnal preferred points. Te main contributions of the introduced algorithm are listed below: (i) We introduce a new RP-based algorithm called RS-KKTPM, based on the S-KKPM metric, by integrating S-KKPM with NSGA-II to obtain the PO solutions in ROI. (ii) Obtaining diferent ranges of ROI in a single run. (iii) Adding fexibility for several RPs at the same time.
(iv) Obtaining excellent performance when the RP is located in diferent regions. (v) Obtaining a good balance between convergence and diversity aspects around the RP. (vi) Solving diferent shapes of PF (e.g., convex, concave, concave, and discontinuous) with a diferent number of objective functions (up to 10 objectives). (vii) Making the results competitive compared with those of the other preference mechanisms on many-objective problems.
Te layout of this work is as follows. Section 2 reviews some fundamental defnitions. An overview of the works relevant to this paper is mentioned in Section 3. In Section 4, the R-NSGA-II algorithm is combined with the S-KKTPM metric. In the following section, the obtained experiments and the results are discussed and described. Section 6 summarizes the paper's achievements and presents some upcoming works. Table 1 displays the nomenclature/abbreviations used in this study.
In an MOP with contradictory objectives, the search space is only partially ordered, and two solutions may be indiferent to each other. A single decision variable simultaneously optimizing all the objectives is unusual. Consequently, for MOPs, the 〈, 〉, and � operators are extended as follows.
Defnition 1 (Pareto dominance relation). Given two solutions x, y ∈ R n , x is said to dominate y in the Pareto sense (denoted by Defnition 2 (non-dominated solution). A solution x ∈ Ω⊆R n (Ω is the feasible space) is said to be non-dominated if and only if there does not exist another solution y ∈ Ω such that f(x)≺f(y).
Te set of solutions in the search space is called the Pareto solution set (PS). In contrast, the set of all non-dominated vectors in the objective space corresponding to the PS is called the PF [38].
Defnition 4 (PS and PF). Te PS is defned as follows: Te corresponding PF is defned as follows: Defnition 6 (ROI). Te ROI is the projection of the set of preferred efcient solutions in the objective space, i.e., is the closest to the RP f(x RP ). δ denotes the radius of the ROI, which is determined by DM.

Related Works
In Section 3.1, KKT optimality conditions are briefy reviewed. Section 3.2 presents R-NSGA-II in detail.

KKT Conditions.
KKT conditions play an important role in optimization theory. Trough these conditions, it is possible to know if the solution produced by the EMO algorithm is the PO solution or not. For the MOP with inequality constraints, KKT conditions are defned as follows [39]: Te parameters β i and c j are called the Lagrange multipliers for the ith objective function and jth inequality  (4) and (6) are called the equilibrium and complimentary slackness equations, respectively. Te conditions stated in equation (5) ensure feasibility for x while the conditions stated in equation (7) ensure that the parameters c j are non-negative. Te conditions stated in equation (8) also ensure that the parameters β i are non-negative, but at least one of them must be non-zero. In the following section, we briefy discuss R-NSGA-II algorithm.

R-NSGA-II Algorithm.
As mentioned in Section 1, classical EMO algorithms mainly aim to develop a fnite number of random solutions into a set of non-dominant solutions that converge and distribute across the entire PF over several generations. On the contrary, preference-based algorithms aim to produce non-dominated solutions centered around the desired part (s) of the PF based on the preference information supplied by the DM. Tis information can be given in several techniques: RPs, aspiration levels, weights, and reference direction [2]. RPs are one of the most used techniques in preference-based EMO algorithms. Usually, an RP is said to be achievable if it lies in the feasible objective space; otherwise, it is said to be unachievable. In 2006, Deb and Sundar [26] put forward R-NSGA-II method, which presents the DM's preferences as one or more RPs. Te method is based on the benchmark manner, which is based on preference information [40]. It is a modifcation of the widely used EMO approach NSGA-II, in which an Euclidean distance metric is applied instead of the crowding distance metric from the RP that indicates DM's preference. Te primary notion behind R-NSGA-II is to give preference to parents who have short Euclidean distances to the RP. Te following is the description of the R-NSGA-II procedure: P t (of size N) is a randomly generated parent population. A new ofspring population Q t (of size N) is generated using the number of operations (binary tournament selection, recombination, and mutation). Tereafter, the populations P t and Q t are combined, and the resulting population R t � P t + Q t (size 2N) is classifed according to dominance in fronts. Te new population P t+1 is built starting with the fronts with the lowest rank until reaching a front F i , which cannot be accepted without making the size of population to exceed N. Next, the preference operator is applied to the front F i to maintain the size of the new population. Te fnal front F i , which cannot be fully accepted, is then considered, and the remaining slots are flled according to an environmental selection approach. Te Euclidean distance for each RP is calculated with respect to each solution of the front F i . For each RP, the solution closest to the said point takes the preferred distance value of 1. Te solutions that are closest to all of the RPs are given the shortest preferred distance. Te preferred distance value of 2 is then applied to the solutions with the next smallest distance to each RP, and the process is repeated for the remaining F i solutions. In the generation of the new population of descendants, the preferred solutions in the selection by the tournament are those with a lower value of preferred distance.
Te idea of ϵ-based selection strategy is utilized to maintain diversity in the solutions close to each RP. First, a solution of the front F i is randomly chosen. Next, the Euclidean distance in the objective space of all the solutions is computed with respect to the chosen solution. After that, the points that have a sum of the normalized diference in the objective search space values less than or equal to ϵ from the selected point are given an artifcial large distance to remove them from the competition; in this method, only a solution within ϵ-neighborhood is relevant. Te process continues by randomly choosing another solution diferent from the previous one, to which the concept of ϵ−based selection strategy described above is applied again.

Advantages and Disadvantages.
Compared with classical RP-based algorithms, R-NSGA-II works well for high-dimensional MOPs; it is suitable for any frontier shape, several objectives, and variables. It also shows some advantages: the classical methods depend on the reference direction (weight vector); however, R-NSGA-II is independent of the weight vector. Moreover, the classical methods in most cases can only fnd efcient solutions for diferent RPs by applying the algorithm to each RP for several times, whereas R-NSGA-II can produce a set of efcient points for diferent RPs in a single simulation. RPs can exist anywhere in the objective space (achievable or unachievable). However, it requires a parameter ϵ to maintain a diversity of selected solutions near the RPs.
As mentioned above, the crowding distance metric of NSGA-II has been replaced using the Euclidean distance metric in R-NSGA-II to obtain the solutions closest to the RPs assigned by DMs. However, the fnal minimal Euclidean distance value is unbeknown. Tus, ascertaining whether the efcient solution is accomplished at the end of an optimization run is difcult. In other words, the Euclidean distance metric does not have any information about the proximity of a solution to the PF. Additionally, in the case of achievable RPs, the Euclidean distance metric may not be monotonically reduced to its minimum value. One major disadvantage of this method is that the DM cannot control the size of each preferred region separately. Furthermore, the DM cannot smoothly control the obtained PO solutions within each desired region. Below, we introduce a new approach that is based on integrating the S-KKTPM metric with the R-NSGA-II algorithm.

The Introduced R-NSGA-II with S-KKTPM
In Section 4.1, the development of the KKT-proximity measure is introduced. Section 4.2 presents the proposed RS-KKTPM in detail.

S-KKTPM.
KKT conditions are necessary to know whether the solution obtained by the EMO algorithm is a KKT point. Hence, they play an important role in optimization theory [39,41]. During the last decade, a KKTproximity measure has been developed utilizing KKT optimality theory. In 2013, a KKT-based proximity metric Computational Intelligence and Neuroscience 5 (KKTPM) was suggested by Dutta et al. [42] to calculate a KKTPM value for any iteration (or solution) x k for a singleobjective optimization problem. Deb and Abouhawwash et al. [37,43] extended the above KKTPM for MOPs. Teir expansion, which is based on the incorporation of the KKTPM metric via scalarization approaches, aims to relate the convergence property of a solution from a specifc optimal solution. Other information on KKTPM for MOPs can be found in [37,43]. In 2021, Eichfelder and Warnow [44] proposed a new KKTPM metric for MOPs without using any scalarization approach. Te authors defned the following methodology for calculating the KKTPM value for any solution x k ∈ R n , for the MOP mentioned in equation (1): where J and M, respectively, are the numbers of constraint and objective functions. Te value ξ k obtained after the optimization is the KKTPM at the point x k . First-order derivatives of constraint and objective functions are necessary to solve this problem. KKTPM metric can be utilized to single, multi, and many-objective optimization problems. Te above problem has (M + J + 1) variables, (M + 2J + 2) inequality constraints, and one equality constraint. To reduce the number of constraints in the optimization problem mentioned above, we propose to redefne it as follows: where the variable vector of the above optimization problem is (ξ k , β, c). Te value of ξ k , which solves the above problem, is referred to as the simplifed KKTPM (S-KKTPM). Te primary goal of reducing the number of constraints is to save the computational cost of an optimization problem. Te above problem has (M + J + 1) variables, (M + J + 2) inequality constraints, and one equality constraint. Te number J of inequality constraints has been reduced compared to the optimization problem mentioned in equation (9) without afecting the optimization process. To ensure that values of both ξ k and ξ k obtained after the optimization are identical at point x k , frst we consider the ZDT1 unconstrained problem with thirty variables [45]. We ran NSGA-II for 200 generations in this problem, with a population size of 40. Computational Intelligence and Neuroscience and ξ k are identical. Second, we consider the SRN unconstrained problem with two variables and two constraints [46]. In this problem, we also ran NSGA-II until generation 500, with a population size of 200.   Computational Intelligence and Neuroscience An advantage of S-KKTPM metric given in equation (10) is that it predicts the convergence behavior of a point from the PF without prior knowledge of the PO solution. Now, we describe several features of the S-KKTPM [37,43,44]: (i) It can be utilized as a termination condition for the algorithm of optimization. (ii) It is applicable in high-dimensional MOPs; S-KKTPM is suitable for any frontier shape, large number of objectives, and variables. (iii) It provides a monotonous characteristic of the S-KKTPM surface over the objective space. S-KKTPM value decreases monotonously almost to zero as the iterate approaches the efcient solution. In this study, we use S-KKTPM optimization problem to calculate ξ k value at iterate x k . We used MATLAB fmincon() algorithm optimization to solve S-KKTPM optimization problem (see Algorithm 1).

Te Proposed RS-KKTPM.
To make R-NSGA-II solutions preferred and acceptable to DMs and to easily control the size of each region, S-KKTPM metric is integrated with the-NSGA-II algorithm.
In this study, we refer to the RP-based S-KKTPM as RS-KKTPM. Te introduced algorithm allows DMs to apply any number of RPs. RS-KKTPM also allows the DMs to control the size of the preferred parts separately. In the introduced algorithm, we replace the Euclidean distance metric, utilized in R-NSGA-II, with S-KKTPM metric. Solutions with small S-KKTPM values are chosen in the introduced method. Te preference operator is utilized in this algorithm to select a subset of solutions from the fnal front that cannot be accommodated totally to maintain the size of population in the novel population. Instead of using the preference distance as in R-NSGA-II, this preference operator uses the preference S-KKTPM metric.
We now characterize an iteration of the introduced R-NSGA-II with S-KKTPM process in which the DM provides one or more RPs in the following section (see Algorithm 2). Both parents and children are merged as usual, and the non-dominated sorting strategy is employed to classify the merged population into non-domination levels (so-called fronts).
Te following are the primary ideas underlying selecting the preferred set of solutions within the preferred range: (i) Solutions closest to RP are always prioritized.
(ii) Preferred-region sizing strategy is used to control the preferred range near RP. (iii) ϵ-based selection strategy is utilized to keep the spread of solutions within the range assigned by the DMs.
Te following changes are made to the original NSGA-II niching approach to integrate the three notions mentioned above: Step 1. Generating a desired region for each RP. Te Euclidean distance between all members from the merged population and an RP is computed to specify the desired region. Ten, the member that has the least Euclidean distance to RP is identifed. Te specifed member (or point) is called mid-point as illustrated in Figure 7.
Step 2. Determining the size of the desired region for each RP. Here, we introduce a new strategy to determine the size of each desired area as follows. Te solution within δ distance of the mid-point is chosen to be in the desired area. Parameter δ is given by the DM, which determines the size of the ROI, as illustrated in Figure 7. Tis fgure also shows how to choose a population of size eight from the merged population containing 17 members. All solutions in the frst front are selected, as shown in Figure 7. Ten, we need only two solutions from the second front. Te remaining two solutions are chosen (from the second front) as follows. Te S-KKTPM value is calculated for each solution (x) within the ROI. Ten, the minimum of the appointed ranks is appointed as the S-KKTPM value to a solution (x). If the solution (x) is not within the preferred region, we set a high value for S-KKTPM (see Algorithm 3). In this manner, the smallest S-KKTPM Input: Population size (N), set of reference points RP← RP 1 , RP 2 , . . . , RP L , Generation number, P Crossover , ϵ Parameter Output: Children (1) Create initial parent population P t of size N; (2) Repeat (3) Generate ofspring population Q t from P t by applying selection, crossover, and mutation operators; (4) Combine P t and Q t population (i.e., R t � P t + Q t ); (5) Classify R t into diferent fronts (F 1 , F 2 , etc., where F 1 is the best non-dominated front, F 2 is the next best non-dominated front, and so on) utilizing non-dominated sorting algorithm; for j← 1 to Ndo (4) dis(RP i , x j ) � Euclidean distance between RP i and x j ; . , x min dis L ; (9) for i← 1 to Ldo (10) for j← 1 to Ndo (11) di s(MP i , x j ) � Euclidean distance between MP i and x j ; (12) end for (13) end for (14) for i← 1 to Ldo (15) for j← 1 to Ndo Calculate S-KKTPM value at iterate x j //Algorithm 1; (18) else (19) Set S-KKTPM (x j ) equal to 10 6 ; (20) end if (21) end for (22)  Computational Intelligence and Neuroscience 9 Step 3. Good distribution of the obtained solutions. Te ϵ-clearing selection strategy, employed in the original R-NSGA-II, is used in RS-KKTPM to control the diversity of chosen solutions near the RPs. A solution is selected randomly from the set of non-dominated solutions to implement this strategy. Ten, any solution with a sum of normalized diferences in objective values less than ϵ is selected and then given a high-preference distance value to discourage it from remaining in the next generations of the evolution process. Te way is then repeated with a new solution picked from the set of efcient points (excluding the one previously selected). Te value of ϵ is selected according to the application and can be diferent for each objective. Tus, it is formed as a parameter provided by the DMs. Figure 8 depicts how to determine the size of the ROI for each RP using the mid-point strategy. As discussed in step 1 above, the mid-point is a member of the population that is closest to the RP. As shown in Figure 8, RP can exist anywhere in the objective region (feasible or unfeasible), whereas the mid-point can exist anywhere in the feasible objective domain only. Te purpose of the proposal of midpoint strategy can be summarized as follows: (1) getting PO solutions that are close to the given RP; (2) determining the size of the ROI by calculating the Euclidean distance between each solution and the mid-point (each distance value is normalized using zero as the lower bound and one as the upper bound to stay within the interval [0, 1]; the solutions that lie within δ value are candidates to be within the ROI); and (3) obtaining a good convergence of solutions towards the ROI. As discussed in step 2 above, the S-KKTPM metric acts as a diferentiator in selecting a solution that should remain in the next generations of the optimization process. Te solution with the smallest S-KKTPM value is preferentially kept for the next generations because it is the closest to the true PF. Tis way, the RS-KKTPM can obtain good convergence of solutions towards the ROI. Te introduced algorithm can well distribute solutions along the preferred part. RS-KKTPM works well with diferent RPs (feasible or infeasible) in the objective space, as displayed in Figure 8. In real-world applications, objectives should be normalized when they do not have the same units. Otherwise, δ is not a meaningful parameter.
One of the essential advantages of the introduced method is its ability to control the size of the preferred areas separately by a single simulation run (see Figure 8). Tis is done using the preferred-region sizing strategy discussed above, based on the S-KKTPM metric. Tis metric is used as a preference operator to select a subset of solutions close to the PF in order to move to the next population. As the iteration approaches the PF, S-KKTPM value decreases monotonically almost to the fnal minimum value (zero). Tis means that the S-KKTPM metric can know the proximity of a point in the search space to the PF. Trough this strategy, the introduced algorithm can steer the solutions during the optimization process towards the preferred regions in proportion to the size of each area. In other words, the large ROI gets more PO solutions compared to the smaller preferred region.
On the other hand, the original R-NSGA-II algorithm cannot control the size of the preferred regions separately through a single run. Te reason is the preferred-region sizing strategy used in this algorithm, which is based on the Euclidean distance metric. Tis metric is utilized as a preference operator in the R-NSGA-II algorithm. However, the Euclidean distance metric does not have the unique properties that the S-KKTPM metric does. For example, the fnal minimal Euclidean distance value is unknown. In other words, the Euclidean distance metric does not have any information about the proximity of a point to the PF. So, the R-NSGA-II algorithm cannot obtain diferent ranges of ROI in a single run.

Experimental Results and Discussion
Tis section uses a set of benchmark problems and engineering design problems to test our introduced methodology. Specifcally, we adopted fve two-objective unconstrained problems taken from the ZDT test suite [45], four bi-objective constrained problems (BNH, SRN, OSY, and TNK) taken from [46], and seven test problems having from three to ten objective functions taken from the DTLZ test suite [47]. In addition, we adopted two engineering design problems, the welded beam design problem with two objective functions (taken from [48]) and the car side impact design problem with three objective functions (taken from [49]). Ten, we compare the performance of the RS-KKTPM approach with six EMO preference approaches, including R-NSGA-II, g-NSGA-II [27], r-NSGAII [28], R-NSGA-III [50], WV-MOEA-P [51], and MOEA/D-PRE.
Te parameters of the suggested method are set as follows: (i) Reproduction operators: as suggested in original study [26], simulated binary crossover (SBX)  Computational Intelligence and Neuroscience probability and SBX index are set to 0.9 and 10, respectively, and the polynomial mutation probability and mutation index are set to 1/n and 20, respectively. (ii) Population size, maximum number of generations, RPs, and size of ROI (δ): diferent parameters for a set of diferent test instances are displayed in Tables 2  and 3.
For constraint handling in constraint test problems and engineering design problems, we handled it by adding a penalty proportional to the constraint violation to the objective function value as suggested in the original NSGA-II algorithm. In minimization problems, this is a popular approach to deal with constraints in evolutionary algorithms.
Te proposed RS-KKTPM algorithm is implemented in the MATLAB R2019a platform. Te source codes for the comparison methods are provided by PlatEMO [51] or downloaded from the authors' home page. Te suggested and compared methods are simulated on a personal computer with an Intel(R)Core(TM)i7-7500 2.9 GHz Quad-Core Processor and 8 GB RAM.

Experiments on Two-Objective Unconstrained ZDT
Problems. Now, we apply our proposed approach to ZDT1 unconstrained problem (it has a convex PF) with thirty variables. Figure 9 illustrates the infuence of diferent values of δ on the distribution of solutions obtained by RS-KKTPM after 200 generations (i.e., 16000 evaluations, given that RS-KKTPM evaluates 80 ofsprings per generation). Tree RPs Trough diferent values of δ, the proposed algorithm can steer the solutions towards the preferred regions in proportion to the size of each region. Parameter ϵ is still required to ensure that the obtained solutions are well distributed within preferred region. In this problem, the parameter ϵ = 0.005 is chosen. Te solutions obtained are clustered near the RPs, as shown in Figures 9(a)-9(c). Te distribution of the obtained PO set depends on the range of each desired region. In particular, the range of solutions obtained is equally vast when the value of δ is large. One of the advantages of RS-KKTPM is that it allows us to adjust the ranges for the desired region in a single run. Tus, if the DM wants to get a set of solutions (near each preferred region) whose number varies depending on the size of each preferred region separately, diferent values of δ can be chosen. In other words, the DM can control the spread of the generated ROIs by changing the value of parameter δ. If δ = 0.5, the RS-KKTPM provides an approximation of the entire PF. On the contrary, Figure 10 shows the PO set produced utilizing R-NSGA-II for the same three RPs on ZDT1 problem. R-NSGA-II is also performed with ϵ = 0.005 and a population of size 80. It is run until 200 generations. Figure 10 shows that the DM (by RNSGA-II) cannot obtain diferent regions of desired regions in a single run. Also, R- NSGA-II cannot steer the solutions toward the preferred regions in proportion to the size of each region. Henceforth, the parameter ϵ � 0.001 is used in all problems. First, we consider ZDT1 test problem with fve RPs, of which three are infeasible and two are feasible, as shown in Figure 11. Each RP and corresponding size of ROI are shown in Table 2. RS-KKTPM is utilized for this problem, where the population members and the maximum number of generations are 40 and 200, respectively. Te parameter is set to 0.05 for each ROI. Figure 11 also demonstrates how easy the proposed algorithm can be modifed to address multiple RPs. As a result, it discovers various ROIs. Well-convergent non-dominated solutions are obtained on PF near all the fve RPs.
ZDT2 is the next problem which has a non-convex PF. Two RPs are chosen, of which one is feasible and the other is infeasible, as presented in Table 2. Te range of each region corresponding to an RP is also presented in Table 2. Te population members and maximum number of generations, respectively, are 40 and 200. Figure 12 displays the convergence and distribution of the solutions near the two chosen RPs. As shown in the fgure, RS-KKTPM algorithm can easily deal with feasible and infeasible RPs. Te proposed algorithm proves its ability to converge and distribute the solutions obtained within the desired ranges provided by the DMs, as illustrated in Figure 12. RS-KKTPM also showed good distribution on this problem when RP is in the infeasible region. 14 Computational Intelligence and Neuroscience Te test problem ZDT3, with 30 variables, has a disconnected set of PFs. Tree RPs are selected (see Table 2), of which one is infeasible and two are feasible. Te desired solutions produced by RS-KKTPM and R-NSGA-II are illustrated in Figures 13 and 14. Te population members were 40, and the maximum number of generations was 200. Tese two fgures demonstrate that our approach is able to steer solutions towards the PF in proportion to the size of each ROI, while the R-NSGA-II cannot. As illustrated in Figure 13, our approach does not get stuck in any locally PO part, and all generated solutions are non-nominated and global PO solutions.
Next, the test problem of ZDT4 with 10 variables is solved utilizing RS-KKTPM and R-NSGA-II. Tis problem has many local PFs. One RP is used with a range of ROI of 0.15, as displayed in Table 2. Te RP is (0.6, 0.6), and the generations are 500. Te plot of the desired solutions produced by RS-KKTPM and R-NSGA-II is represented in Figures 15 and 16, respectively. As illustrated by the two fgures, the performance of RS-KKTPM is much better than that of R-NSGA-II in terms of the distribution and convergence of solutions from the PF. As shown in Figure 15, the selected RP is somewhat far from the PF, which indicates the ability of the introduced approach to work well in the case of distant RPs. Tus, as shown in Figure 15, even though the problem has more than 100 local fronts, the introduced algorithm can converge well to the true PF.
Finally, we apply our proposed method to a ZDT6 problem that has a non-convex PF. Figures 17 and 18  respectively. Both techniques used the same RPs, the same number of population members, and the same number of generations (see Table 2). Tree RPs are chosen, of which the frst lies in feasible search space, the second lies close to/on PF, and the third lies in infeasible search space.  Tis means that if the DM wants to get PO solutions of diferent sizes for all regions, the introduced algorithm can do that. In contrast, R-NSGA-II cannot control the number of solutions for each desired area. Tis is because parameter ϵ takes only one value for all preferred regions corresponding to the given RPs. In other words, when multiple RPs exist, R-NSGA-II cannot give diferent values for ϵ in a single run. Tis means that if the DM wants to get PO solutions of diferent sizes for all regions, the R-NSGA-II algorithm cannot do that.

Experiments on Two-Objective Constraint Problems.
We now consider two-objective constraint problems: BNH, SRN, OSY, and TNK [46]. RPs and some essential parameters used to solve these problems are shown in Table 2. BNH, TNK, and SRN have only two constraints and two variables. First, the efcient solutions obtained by the RS-KKTPM on BNH with three RPs are illustrated in Figure 19. It is clear from this fgure that our approach can fnd the desired regions near the RPs. Second, the solutions obtained on SRN with two RPs are displayed in Figure 20. Te RS-KKTPM algorithm works well when the RP is in the feasible or infeasible domain, as displayed in Figure 20. Next, we consider the OSY test problem, which has six constraints and six variables. Figures 21 and 22 show the obtained solutions by RS-KKTPM and R-NSGA-II on OSY, respectively. Two RPs are chosen with a range of ROIs and a population � 40 (see Table 2). Although RS-KKTPM cannot converge to the true PF, it converges slightly better than R-NSGA-II, as shown in Figures 21 and 22.
Finally, the desired regions obtained by the RS-KKTPM and R-NSGA-II on the TNK problem are illustrated in Figures 23 and 24, respectively. Two RPs and population members are chosen as displayed in Table 2. As shown in  Figures 23 and 24, the performance of our approach is a little similar to R-NSGA-II. In summary, the introduced approach balances diversity and convergence around ROI for constraint test problems and handles any number of predefned RPs.

Experiments on Tree-Objective Problems.
We will select the original DTLZ1, DTLZ2, and DTLZ5 and their scaled versions. Table 3 provides some information about these problems and some parameters required. First, the DTLZ1 problem contains many local PFs, possibly causing some points to stop. Tis scenario is a relatively complicated problem to address for global optimality. Figure 25 shows the obtained preferred PO solutions using RS-KKTPM and R-NSGA-II algorithms on the three-objective DTLZ1 (DTLZ1 3 ) problem, and the parameter values are presented in Table 3. Te three aspiration points are chosen. Te distribution and convergence of solutions found by RS-KKTPM are substantially superior to those by R-NSGA-II, as displayed in Figure 25. Next, RS-KKTPM is utilized to solve the three-objective DTLZ2 problem. Figure 26 shows the obtained solutions with 60 population members with two RPs. Figure 26 clearly illustrates that the RS-KKTPM algorithm can access the efcient region of true PF with few number of population sizes, thereby helping the DM determine the required ROI easily. Finally, the RS-KKTPM is utilized to solve the three-objective DTLZ5 problem. Te two RPs, RP 1 � (0.6, 0.6, 0.65) T and RP 2 � (0.2, 0.3, 0.8) T , are used. Te sizes of the preferred areas for RP 1 and RP 2 are 0.2 and 0.1, respectively, as shown in Table 3. Our algorithm is employed to solve this test problem with 60 populations and runs up to 300 generations. Te obtained preferred areas of the true PO solutions are displayed in Figure 27. Te obtained solutions are distributed according to the size of each preferred area. In a single simulation run, both areas are discovered. Note that the number of solutions generated in the frst preferred area, corresponding to RP 1 , is greater than that generated in the second preferred area, corresponding to RP 2 . Well-convergent and well-distributed solutions are obtained on PF in all two RPs. Tus, the DM can control the size of each efcient region (s) of the true PF separately and in a single simulation run.

Experiments on Many-Objective Problems.
Finally, we test our introduced approach on the many-objective versions of the problems of DTLZ1 and DTLZ2. Table 3 displays all parameters used for 5 and 10-objective problems. First, RS-KKTPM is used for DTLZ1 5 and DTLZ1 10 problems. Population sizes of 80 are used for the two problems. Figures 28 and 29 present the obtained part in a parallel coordinate plot. One RP is used for each problem, as shown in Table 3. Te PO solutions of these problems must satisfy M i�1 f i � 0.50. RS-KKTPM can discover the needful regions of the efcient set corresponding to the one predefned RP by the DM.
Finally, the RS-KKTPM algorithm is applied for DTLZ2 5 and DTLZ2 10 with 14 and 19 decision variables. Tis algorithm is applied for these test problems with 60 populations and runs up to 300 generations. Figures 30 and 31 show the obtained solutions for DTLZ2 5 and DTLZ2 10 with two and one aspiration points, respectively. Te PO solutions to these problems must obey the next equation

18
Computational Intelligence and Neuroscience

Welded Beam Design
Problem. First, we now employ a two-objective welded beam design problem [48] as a realworld example. Te frst objective is to minimize the cost of fabrication, whereas the other objective is to minimize the end defection of the welded beam. Te design of welded beam structure is shown in Figure 32. Tis problem involves four decision variables, namely, h (weld thickness), l (clip length), t (the height of bar), and b (the thickness of bar). It has also four non-linear constraints. Te problem is mathematically formulated as follows [48]: minimize: where υ � 6 * 10 3 (sqrt(2) * hl) ,  Figure 33(a) shows that the introduced algorithm outperforms the others by adjusting the size of each preferred region (separately) corresponding to the supplied RP. It has the ability to steer the solutions towards the preferred regions in proportion to the size of each area. As shown in Figure 33(a), the obtained solutions in the preferred region, corresponding to RP 1 , are more compared to the obtained solutions in the other two preferred regions. Indeed, the size of the ROI, corresponding to RP 1 , is greater than the sizes of the other two preferred regions, i.e., δ 1 > δ 2 and δ 1 > δ 3 . On the other hand, the preferred region, corresponding to RP 3 , contains a few solutions compared to the other two preferred regions because δ 3 < δ 1 > and δ 3 > δ 2 . Figure 33

Car Side Impact Design Problem.
Car side impact design is a constrained optimization problem [49]. Tis problem has three optimization objectives which are described as follows: the frst is to reduce the car's weight, the second objective is to minimize the pubic force experienced by a passenger, and the last objective is to minimize the average velocity of the V-Pillar responsible for withstanding the impact load. It has seven decision variables: B-Pillar, door beam, B-Pillar inner reinforcement, foor side inner, door beltline reinforcement, cross members, and roof rail (see Figure 34). Te mathematical model of this problem is as follows [49]: minimize: f 1 (x) � 1.98 + 4.9v 1 + 6.67v 2 + 6.98v 3 + 4.01v 4 + 1.78v 5 + 0.00001v 6 + 2.73v 7 , (13) Figure 35 shows the solutions produced by RS-KKTPM, R-NSGA-II, R-NSGA-III, and MOEAD-PRE on the car side impact design problem. In this problem, all relevant parameters of the four comparison approaches are briefy presented below. Two RPs are chosen: RP1 � (40, 3.5, 11) and RP2 � (26, 4, 11.5). For this test problem, the comparison algorithms are used with 80 population members and run until 300 generations. For the RS-KKTPM algorithm, the radius corresponding to each RP is δ 1 � 0.1 and δ 2 � 0.05. For the rest of the algorithms, the size of preferred regions is equal to 0.1. In this problem, all objective function values are normalized using the ideal point (15,3,10) and nadir point (50,5,14). Figure 35(a) shows that the suggested algorithm can control the number of solutions for each desired area in proportion to its size. In other words, a larger desired area gets more solutions, while a smaller preferred region gets fewer solutions. As displayed in Figure 35(a), the number of solutions for the frst preferred area, corresponding to the P 1 , is more than that for the second desired area, corresponding to the RP 2 . Tis is because the size of the frst region is greater than the size of the second region, i.e., δ 1 > δ 2 . Terefore, RS-KKTPM can control the size of each desired area separately, while the rest of its algorithms cannot do that (see Figures 35(b)-35(d)). Tus, if the DM is interested in fnding solutions in region of diferent sizes, the introduced algorithm can fnd solutions near the RPs and proportional to the size of each preferred area separately.

Performance Metrics.
No single performance metric can provide an accurate assessment of an EMO's performance [54]. In our empirical investigations, we use two of the most recognized performance metrics to assess the quality of preferable efcient solutions of preference-based EMO algorithms: R-HV and R-IGD. [55]. Both metrics are utilized to detect the ROI's convergence and the diversity of efcient solutions simultaneously. Tey are based on two performance metrics, the hypervolume (HV) metric and the inverted generational distance (IGD) metric, which are designed for the entire PF and applicable for partial preferable efcient solutions. Te larger the R-HV values are or   As mentioned earlier, S-KKTPM requires the gradient of all objective and constraint functions. Ten, algebraic calculations are performed to compute the theoretical closeness of x to the true optimal solution. For MOPs, S-KKTPM calculates the closeness metric from a specifc PO point. In this article, the introduced RS-KKTPM approach is then compared with six EMO approaches. For a fair comparison between all algorithms, we used equal function evaluations in the comparison.
Compared to the number of function evaluations required for a solution evaluation, the savings reported for the S-KKTPM calculation may not be signifcant because the evaluation of S-KKTPM for all solutions is an additional computational expense and requires more computation. Tus, once gradients are computed for a real-world problem, the computational time needed for the S-KKTPM optimization procedure would make a small addition to the overall computational time. In the meantime, S-KKTPM helps improve convergence and can diferentiate between diferent non-dominated solutions that are not applicable by using the Euclidean distance or any other evolutionary algorithm.
Tables 5 and 6 display the mean and standard deviation of R-HV and R-IGD values, respectively. Te best mean of R-HV and R-IGD metrics is highlighted in bold in Tables 5 and 6.
According to R-HV metric, the ROI is approximated by the introduced RS-KKTPM algorithm in a better way than other algorithms for all the examined problems except the ZDT4, DTLZ2 3 , and DTLZ2 5 test cases (see Table 5). We obtain almost the same results according to the R-IGD metric, as shown in Table 6. Based on the R-HV values and the R-IGD values, RS-KKTPM demonstrates better distribution and convergence than other algorithms. Te practical fndings on the 14 benchmark test problems illustrate that the RS-KKTPM approach outperforms the other approaches used in 11 of 14 comparisons.

Conclusions
In this study, the RS-KKTPM preference-based EMO algorithm is proposed. It is an expansion of the R-NSGA-II method, where the Euclidean distance metric is replaced by the S-KKTPM metric. Te following are the properties of this new algorithm: (i) Te RS-KKTPM can obtain the ROI at any specifc position of the RP (in the feasible area, on/near the PF, and the infeasible area). (ii) Te range of each obtained ROI can be controlled by adjusting the interest radius size of each ROI separately and in a single simulation run.  24 Computational Intelligence and Neuroscience (iii) Te RS-KKTPM algorithm, given herein, improves the quality of the PF approximation and allows a uniform distribution of the approximating objective vectors. (iv) Te performance of RS-KKTPM is better than that of R-NSGA-II, g-NSGA-II, r-NSGAI, R-NSGA-III, WV-MOEA-P, and MOEA/D-PRE on most multi and many-objective problems.
Te direction of future research focuses on using the S-KKTPM metric to improve the performance of other EMO optimization algorithms by reference direction approaches, such as MOEA/D and NSGA-III. Tese approaches can also be utilized to solve engineering design problems and highly complex problems.

Data Availability
No data were used to support this study.

Conflicts of Interest
Te authors declare that they have no conficts of interest.