An enhanced surrogate-assisted differential evolution for constrained optimization problems

The application of evolutionary algorithms (EAs) to complex engineering optimization problems may present difficulties as they require many evaluations of the objective functions by computationally expensive simulation procedures. To deal with this issue, surrogate models have been employed to replace those expensive simulations. In this work, a surrogate-assisted evolutionary optimization procedure is proposed. The procedure combines the differential evolution method with a k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k$$\end{document}-nearest neighbors (k-NN) similarity-based surrogate model. In this approach, the database that stores the solutions evaluated by the exact model, which are used to approximate new solutions, is managed according to a merit scheme. Constraints are handled by a rank-based technique that builds multiple separate queues based on the values of the objective function and the violation of each constraint. Also, to avoid premature convergence of the method, a strategy that triggers a random reinitialization of the population is considered. The performance of the proposed method is assessed by numerical experiments using 24 constrained benchmark functions and 5 mechanical engineering problems. The results show that the method achieves optimal solutions with a remarkably reduction in the number of function evaluations compared to the literature.


Introduction
EAs have become appealing to solve complex real-world engineering optimization problems, for which classical mathematical methods may no longer be effective. This includes, for instance, the design of offshore oil and gas production systems such as floating platforms, risers, mooring systems and submarine pipelines (Vieira et al. 2012;Lucena et al. 2014;Monteiro et al. 2016): those problems are characterized by many design variables, defined on high-dimensional spaces, with objectives and constraints that are generally nonlinear functions. However, the use of EAs for such problems can become impractical due to excessive computational costs: many evaluations of candidate solutions may be needed until a satisfactory solution is found, and each evaluation requires computationally expensive numerical simulations.
Those issues can be handled by building cheaper approximation models (usually referred as surrogate or metamodels) (Queipo et al. 2005) to partially or totally replace the numerical simulation procedures, emulating or mimicking the behavior of the actual model as closely as possible by capturing the relationship between a given set of inputs and outputs. The goal is to significantly reduce computational costs, while presenting accurate results. In this sense, several surrogate models have been proposed. They may comprise polynomial approximation models (Saunders et al. 1998;Hendrickx and Gorissen 2006); radial basis functions (Kybic et al. 2002a(Kybic et al. , 2002bMullur and Messac 2006;Lowe 1988); response surface models (Myers et al. 2016); artificial neural networks (Ferrari and Stengel 2005); Gaussian process (Emmerich et al. 2006;Ulmer and Streichert 2003); Kriging (Simpson et al. 2001); support vector machines (Kecman 2005); and similaritybased models such as the nearest-neighbor (NN) algorithms (Krempser et al. 2012(Krempser et al. , 2017Fonseca et al. 2009). Particularly, the k-NN model presents many advantages in terms of good performance and ease of implementation, as reported in Fonseca et al. (2009); Liu and Sun 2011;Miranda-Varela and Mezura-Montes 2016).
Many applications of surrogate models with EAs have been presented in the literature. In (Grefenstette and Fitzpatrick 1985) a polynomial model was used to approximate the objective function. In (Ratle and Optimal sampling strategies for learning a fitness model. 1999) a kriging metamodel, trained using few samples evaluated by the actual objective function, has been successfully used to reduce the number of function evaluations. The goal of reducing computational costs was also sought in Hong et al. (2003) where the fitness function has been approximated by an online multi-layer neural network, and in Sun et al. (2014) that presented a surrogate-assisted PSO algorithm. Surrogate models have also been employed in association with the genetic crossover operator (Anderson and Hsu 1999), to approximate the objective function (Rasheed and Vattam 2002), and to calculate approximations for the sum of violation terms (Miranda-Varela and Mezura-Montes 2016). The IEEE Congress on Evolutionary Computation (CEC) competitions have presented many articles that solve computationally expensive functions by employing metamodels with EAs. Regarding particularly the application to offshore oil & gas production systems mentioned above, surrogate models based on neural networks were used in Pina et al. 2014a;Pina et al. 2014b) to avoid expensive nonlinear dynamic Finite Element analyses.
Among evolutionary algorithms, the differential evolution (DE) algorithm (Storn and Price 1997) has been widely used due to its outstanding performance. The association of DE with approximation models has been studied in some works, including neural networks (Wang et al. 2010) and neural networks of radial basis functions (Pahner and Hameyer 2000). In (Mallipeddi and Lee 2015) a simple kriging metamodel was employed to comprise an ''evolving surrogate model-based differential evolution'' (ESMDE). As the population evolves over generations, the surrogate model also evolves using the appropriate parameter setting during different stages of the evolution, thus assisting the DE in generating competitive offspring. A surrogate model for DE consisting of multiple local surrogate models was presented in Jin et al. (2015), built from the feedback of the evolutionary process in diverse overlapped local regions of the search space. In (Liu et al. 2016) a multi-fidelity surrogate model-based DE optimization strategy was developed, coupled to a data mining technique to handle the discrepancy between the low and high fidelity simulation models. A surrogate model to assist the DE algorithm in generating competitive solutions along the evolutionary process was proposed in Awad et al. (2018), using an adaptation scheme to adapt parameters of a Kriging model. Recently, a two-layer adaptive surrogateassisted DE has been proposed in Yang et al. (2019), where three different search strategies are adapted during the evolutionary process according to the feedback information, to measure the status of the algorithm approaching the optimal value. Local similarity-based surrogate models based on an r-nearest neighbors have been employed in Krempser et al. (2012) and (Krempser et al. 2017) to solve complex structural optimization problems; the results indicated that the performance of DE has been significantly improved. In this context, in Garcia et al. 2017a we had presented preliminary studies leading to the optimization procedure called SADE-kNN (surrogate-assisted differential evolution with the k-NN surrogate model). Now, this work presents and evaluates an improved surrogate-assisted evolutionary method, intended to comprise a more efficient optimization tool to solve complex real-world engineering problems. The proposed method is comprised by two variants, referred to as SADE-kNN-MCR and SADE-kNN-MCR 1 . Both incorporate a merit-based scheme to manage the database used by the k-NN surrogate model to store the solutions evaluated by the exact model, privileging solutions able to better approximate new candidate solutions. Compared to the authors' previous work (Garcia et al. 2017a), the proposed method introduces novel characteristics, beginning with its association with the multiple constraint ranking (MCR) constraint-handling technique (Garcia et al. 2017b). The MCR comprises a rank-based approach that builds multiple separate queues, corresponding to the values of the objective function and the violation of each constraint. It has been devised to deal with complex engineering problems characterized by a large number of constraints that may present different units and/or different orders of magnitude, where the possible range of minimum and maximum violation values may not be known a-priori (thus standard normalization procedures cannot be applied). Moreover, results of preliminary studies indicated that the previous SADE-kNN method presented by the authors in Garcia et al. 2017a is susceptible to premature convergence and/or stagnation around local optima. This is a known issue related to DE algorithms, as acknowledged in the literature (Lampinen and Zelinka 2000;Hrstka and Kučerová 2004). To circumvent this issue, the variant called SADE-kNN-MCR 1 incorporates a strategy that triggers a random reinitialization or ''rebirth'' of the population whenever it stagnates on local optima.
The performance of the proposed method is compared with other surrogate-assisted procedures recently presented in the literature, by applying them to several benchmark functions (including the G-suite from the IEEE competition on real parameter constrained optimization (Liang et al. 2006), and mechanical engineering problems), and assessing their relative performance by direct comparisons and nonparametric statistical tests.
This paper is organized as follows: Initially, Sect. 2 summarizes the differential evolution algorithm, the k-NN similarity-based metamodel, and a review of some surrogate-assisted methods for constrained problems presented in the literature. Section 3 describes the variants of the proposed surrogate-assisted differential evolution method (SADE-kNN-MCR and SADE-kNN-MCR 1). The full sets of numerical experiments are described in Sect. 4. The results are presented in Sects. 5 and 6. The former compares both variants of the proposed method with the original SADE-kNN previously presented by the authors in Garcia et al. 2017a. Then, Sect. 4 presents comparisons with other methods recently presented in the literature that are listed in Sect. 0. Final remarks and conclusions are presented in Sect. 7.

Differential evolution
The DE algorithm, proposed by Storn and Price (Storn and Price 1997), has been shown to be a very effective algorithm to solve complex real-world optimization problems such as system design (Storn and System design by constraint adaptation and differential evolution. 1999) and robot manipulator design (Bergamaschi et al. 2008). One of the main characteristics of DE is its simplicity and ease of implementation. A detailed review of its basic concepts is provided by Das and Suganthan (2011);Wu et al. 2018), along with a survey describing its major variants, covering several aspects related to its application to multiobjective, constrained, large-scale and uncertain engineering optimization problems.
Four DE variants were proposed in Kenneth (1999); Price et al. 2006), characterized by the procedures adopted for the selection of candidate solutions, and types of recombination. In this work the first variant (DE/rand/1/ bin) will be adopted. Initially, vectors x i (i = 1,…N) that comprise the individuals of the population are randomly created. Then, the mutation operator defined by the following expression is applied: This operator provides the ''mutant vector'' u i , given by the perturbation of a ''parent'' vector x r 1 with the difference between two other vectors x r 2 and x r 3 (all randomly selected from the population such thatx r1 6 ¼ x r2 6 ¼ x r3 6 ¼ x i ), affected by the parameter F that is a user-defined scale factor.
A crossover operator is also defined, to generate trial vectors v i that combine the solutions in the current population x i and the mutant solutions u i generated by Eq. 2.1. Considering that D is the dimension of each vector, and taking a user-defined crossover rate CR 2 0; 1 ½ ; this operator is expressed as: where d = 1,2,…,D; rand(d) is a random number between 0 and 1; and drand is an integer value randomly selected in the range [1,D] to ensure that at least one variable of the mutant vector u i will be copied into the trial vector v i . After crossover, a greedy selection is performed by comparing the objective function values of solutionsv i .and x i ; only the best one is maintained in the population.

The k-NN similarity-based metamodel
Similarity metamodels evaluate a solution by calculating the average of the objective function of the nearest solutions weighted by its distances. This technique can be classified in two categories, according to the criterion that defines the candidate solutions able for selection: • k-nearest neighbors (k-NN) (Shepard 1968) that consider the k nearest candidate solutions; • r-nearest neighbors (r-NN) (Hu and Lee 2006) that consider a given neighborhood controlled by a range r.
In the k-NN similarity-based metamodel, which is adopted in this paper, the solutions ðx; f ðxÞÞ evaluated by the exact function are stored in a database as an ordered set D. Any solution x generated during the evolution is evaluated by the following weighed average: An enhanced surrogate-assisted differential evolution for constrained optimization problems 6393 where dist x; x a ð Þ is the Euclidean distance between x and x a , x a is a solution stored in the database D, and k is the number of neighbor solutions.
The number of neighbors k and the size of the database D significantly impact the computational costs of the algorithm, since they determine the number of operations to perform the similarity comparisons and to calculate the approximations. Section 1.1 will present some considerations regarding the database management.
As mentioned in Sect. 1, various kinds of surrogate models have been employed in association with EAs. Following the previous SADE-kNN method presented in Garcia et al. 2017a, in this work we will continue using the k-NN approximation to comprise the improved SADE-kNN-MCR and SADE-kNN-MCR 1 methods. This is due to several reasons: initially, because this model has been widely recognized in the literature for its simplicity and good performance (Fonseca et al. 2009;Liu and Sun 2011;Miranda-Varela and Mezura-Montes 2016). Then, because it is particularly suited to the association with the merit-based database management that will be described shortly in Sect. 1.1: while standard k-NN implementations require additional computation costs to assemble the database and perform the neighborhood comparisons, the adopted database management procedure keeps the database size constant, storing only the solutions that provide most accurate approximations. This helps to maintain an adequate balance between accuracy and computational efficiency. Finally, another reason to use the k-NN model is that it is also particularly suited to the ''rebirth'' technique that will be described in Sect. 1.2 to avoid stagnation around local optima. It will be seen that this association enhances the main feature of the k-NN: the use of similarity for approximations, to benefit solutions more similar to the current population, thus increasing the accuracy of the model.

Constraint-handling techniques
A general constrained optimization problem may be formally defined in a D-dimensional search space as follows: Þis the vector of design variables that characterize a given candidate solution, with components x d presenting lower and upper bounds lb d ; ub d ½ . The goal is to minimize the objective function f ðxÞ, considering inequality and equality constraints (respectively g j ðxÞ 0 and h j ðxÞ ¼ 0) that define the feasible region.
Many classical constraint-handling techniques (CHTs) transform this constrained problem into an unconstrained one. For this purpose, whenever any given constraint is violated a penalty function pðxÞ is added to the objective function f ðxÞ; thus, the objective function and constraints are combined into a single value. Another approach is to maintain those values apart along the evolution. Deb (Deb 2000) proposed a binary tournament selection method that ranks the solutions according to the values of the objective function and constraint violation. That method comprises a pairwise comparison where a feasible solution is always preferred over an infeasible one; if both solutions are feasible, the one having the best objective function value is preferred; and if both solutions are infeasible, the one having the lowest constraint violation value is preferred. Later, other methods have followed this rank-based approach, such as the Stochastic Ranking (SR) method (Runarsson and Yao 2000), the GCR (Global Competitive Ranking) method (Runarsson 2002), and the Ho and Shimizu ranking method (HSR) (Ho and Shimizu 2007).

Brief review of some methods
Surrogate-assisted evolutionary methods to solve constrained problems have been presented in the literature, including the ASRES -Approximated Stochastic Ranking Evolution Strategy, which employs two levels of approximation along the evolution, to estimate the values of the objective and penalty functions (Runarsson 2006).
The SADE-kNN method presented in Garcia et al. 2017a employed the binary tournament selection technique proposed by Deb (Deb 2000). Other methods recently presented in the literature also employed the DE algorithm associated with some of the CHTs described above, including for instance SA-DECV (Miranda-Varela and Mezura-Montes 2016) and SA-DECV_SR (Miranda-Varela and Mezura-Montes 2018). These works presented methods that couple DE with the k-NN surrogate model to approximate both the values of the objective function and the sum of constraint violations. The former also employs Deb's feasibility rules of Deb (2000); the latter work (Miranda-Varela and Mezura-Montes 2018), after comparing different CHTs associated to the same surrogate model, concludes that the best performance is provided by the stochastic ranking (SR) method (Runarsson and Yao 2000).
Recently, Yu et al. 2019 presented the e-DE method following two different approaches for the crossover of solutions (with the goal of balancing population diversity and fast convergence). This method employs seven comparison criteria for feasible and unfeasible solutions, comprising an extension of Deb's feasibility rules.

Description of the proposed method
This section describes the main characteristics of the new surrogate-assisted evolutionary method proposed in this work. The method is comprised by two variants, referred to as SADE-kNN-MCR and SADE-kNN-MCR 1 . Those characteristics are related to the following aspects: Both variants include the multiple constraint ranking (MCR) technique, described in Section 3.1, comprising a specialized CHT to deal with the large number of complex constraints that may be found in real-world engineering problems. Both variants also include the merit-based strategy described in Sect. 1.1 to manage the database that stores the exact solutions employed by the k-NN surrogate model. Finally, the variant SADE-kNN-MCR 1 adds the strategy described in Sect. 1.2 that triggers a random reinitialization or ''rebirth'' of the population whenever it stagnates on local optima.

The multiple constraint ranking technique (MCR)
Recently, focusing on the solution of complex real-world engineering optimization problems, we had introduced the multiple constraint ranking technique (MCR) (Garcia et al. 2017b) as a CHT working in association with a genetic algorithm. Now, in this work we will derive an association of MCR with the DE algorithm. As mentioned before, this is one of the main novelties of the variants of the proposed method (SADE-kNN-MCR and SADE-kNN-MCR 1) compared to the original SADE-kNN method described in Sect. 2.3.2 that employed as a CHT the pairwise tournament method of Deb (2000). Section 2.3.1 have described some CHTs that follows a rank-based approach, including HSR (Ho and Shimizu 2007) that employed a single rank R / for all violation values (built as a single sum considering violations in all constraints). MCR extends this approach by splitting R / into several ranks R j / , one for each constraint j ¼ 1; . . .; p, and builds multiple queues considering the objective function values and the violations of each constraint. In MCR, the fitness function F for each solution x is evaluated as follows: where R f ranks the value of the objective function f, and R Nv ranks the number of constraints violated. Algorithm 1 presents the pseudocode of MCR. As seen in row 1, firstly N individuals of the population are randomly created. Then, all individuals have their objective function and constraint violations evaluated (rows 2-4). The fitness function of each individual is evaluated in rows 10 or 12 (depending on whether there are feasible individuals in the population or not) by summing their ranking positions R Nv , R j / and R f that have been calculated in rows 6, 7 and 9 respectively. It is interesting to remark that the two ''for'' loops of rows 2 and 5 cannot be merged into one, since the values of the objective function and constraint violations for the entire population are required before calculating the ranking positions at rows 6, 7 and 9.
As stated in Garcia et al. 2017b, HSR and other previous rank-based techniques may face difficulties in complex problems characterized by constraints with different units and/or different orders of magnitude. This is the case of many practical applications in the offshore oil industry, such as the optimization of risers connected to floating platforms (Vieira et al. 2012;Lima et al. 2005;Pina et al. 2011), or the optimization of subsea pipeline routes (Lucena et al. 2014; Baioco and de Lima Jr. MHA, Albrecht CH, de Lima BSLP, Jacob BP, Rocha DM 2018). For instance, in this latter application many disparate constraints are specified, such as declivity (measured in degrees); minimum radius of curvature (in meters); structural constraints such as on-bottom stability (DNVGL 2017a) and fatigue due to vortex-induced vibrations (VIV) in free spans (DNVGL 2017b). Constraints may be defined even as non-dimensional quantities, such as the number of identified interferences of the route with seabed obstacles (subsea equipment, flowlines, other pre-existent pipelines, regions with corals or geotechnical hazards). In such cases, if just one single rank R / of constraint violation values is built (as in the case of HSR), it could be dominated by a given constraint with higher violation values.
Of course, this wouldn't be an issue for simpler cases where one knows in advance the ranges of minimum and maximum violation values, and normalization procedures could be applied. However, for many real-world engineering applications such as those offshore applications listed above, this information is not available and cannot be estimated in advance. With MCR, all types of constraints are given the same priority and importance in the evaluation of each solution, regardless of the range of their violation values.
MCR presents other advantages: it does not require any user-defined parameter; also, when there are no feasible solutions in the population, F depends only on the ranks associated to the constraints ðR Nv and R j / Þ; thus, the computational costs associated to evaluations of the objective function f are avoided.

Surrogate model: Merit-based database management
A difficult and nontrivial task associated to the k-NN metamodel described in Sect. 2.2. is the definition of an adequate database size D to store the exact solutions. Extensive databases may slow down the evolution; many similar individuals can be generated, compromising the population diversity. On the other hand, small databases may not efficiently approximate the different solutions scattered through the search space. This requires an efficient database management considering the strategies for insertion, maintenance, and replacement of the solutions to keep the database size constant. Usually, an insertion step choses the best solution of the population for evaluation by the exact model (Regis 2014;Regis and Shoemaker 2007;Elsayed et al. 2014). In this work, a merit-based database management strategy is employed to define how the solutions are replaced according to their contribution to the approximation, following some guidelines that have been suggested in Silva (2016). Initially, a predetermined number of solutions are generated and evaluated by the exact function to compose the database D. Those with the best objective function values also compose the initial population N (here it is assumed that N B D). Evolution occurs considering the current population evaluated by the surrogate model. The database solutions used in the approximation process are scored according to their selection as the k-nearest neighbors: the selected solutions have their score increased, while the others have their score decreased.
After each generation, a given number of lower scored solutions (corresponding to a user-defined parameter %N as a percentage of the solutions on the current population) are replaced by solutions evaluated using the exact function and introduced into the database.
Algorithm 2 shows the general structure of the metamodel including this merit-based strategy. The main loop of the algorithm begins at row 1 for all solutions that comprise the population. Then, the loop shown at rows 2-4 calculates the Euclidean distance between the current solution x i and each exact solution x a stored in the database D. The database is then sorted according to those distances (row 5), and the k-nearest neighbors of x i (selected in row 6) are used to calculate the objective function of the current solution x i according to Eq. (2.3) (row 7). Next, the score of all solutions stored in the database D is updated in rows 10 and 13, respectively, to increase the score of the selected k-nearest neighbors, and to decrease the score for the other solutions that have not been selected. Finally, at row 16 after the main loop, a given percentage of the lower scored solutions are replaced in the database for solutions that have been evaluated using the exact function.
This process enhances the main feature of the k-NN model: the use of similarity for approximations. It keeps a memory of the most promising regions, and, while allowing new regions of the search space to be visited (and benefiting from them), there remains a tendency to return to the best regions already known. As a result, this similarity process benefits solutions more similar to the current population, thus increasing the accuracy of the model.

Avoiding stagnation around local optima
A very common problem associated to optimization algorithms in general is the stagnation of solutions around a local minimum in the search space. For a wide class of problems, the distribution of the solutions in the search space may be a decisive factor for this stagnation behavior, when many evaluations of the objective function may be uselessly performed.
An example of this behavior may be observed in Table 1 that presents results of the application of the SADE-kNN-MCR to the G1 function (from the G-Suite benchmark functions proposed for the IEEE-CEC 2006 competition on real parameter constrained optimization (Liang et al. 2006)). The success rate (SuccRate) of the algorithm is 0.84; that is, it finds the optimal solution -15.0000 in 21 out of the 25 runs. In those 21 successful runs, the algorithm performs an average of 3,783 evaluations of the objective function (succ_nfe). This value corresponds to approximately 0.76% of the maximum number of evaluations (500,000) that defines the stopping criterion of the algorithm. In the other 4 runs, the algorithm remains stagnated in the local optimum -12.0000 until it reaches the termination criterion of 500,000 evaluations. This way, the average number of evaluations considering all 25 runs (nfe) increases to 83,177.22. This does not accurately represent the performance of the algorithm: although it was able to reach the optimum with a small number of evaluations in most of the runs (21), the stagnation in the other 4 runs drastically increases the value of this metric (nfe).
There are several strategies to avoid this problem; one of them consists in restarting the population along the evolution, generating new candidate solutions whenever the population converges to a local optimum. The strategy employed in this work maintains in the database the population from the previous generation and triggers the reinitialization (or ''rebirth'') of the population after a given number of generations in which stagnation of the best solution is observed.
In the context of the SADE-kNN-MCR algorithm, this strategy may be particularly suited to avoid frequent premature convergence behavior. Firstly, because the new population may visit new regions of the search space, counteracting the tendency of being stuck in regions where it previously stagnated. Also, because it does not introduce additional computational costs since the objective function of the new individuals from the restarted population are approximated by the k-NN metamodel (using the values from the solutions already stored in the database), therefore not requiring costly extra evaluations of the exact function.

Summary of the proposed method
To summarize the method proposed in this work that combines the DE algorithm with a k-NN surrogate model and the MCR constraint-handling technique, Algorithm 3 presents the pseudocode that describes its variants SADE-kNN-MCR and SADE-kNN-MCR 1 . Both methods begin by randomly initializing (row 1) and evaluating (row An enhanced surrogate-assisted differential evolution for constrained optimization problems 6397 2) the objective function of the individuals stored in the database D. Then, as seen in row 3, the N individuals with the best objective function values are selected to comprise the initial population. The evolutive process is defined by the loop beginning in row 4. At each step, while the maximum number of objective function evaluations is not reached, firstly the DE operators defined in Sect. 0 are applied to generate the new population (row 5). Then, the objective function f ðxÞ and fitness function F ðxÞ are evaluated according to Algorithms 2 and 1 (rows 6 and 7). For the SADE-kNN-MCR 1 variant, the algorithm checks if the population has stagnated around a local optimum (row 9). If this is the case, the current population is randomly reinitialized (row 10) according to the ''rebirth'' procedure described in Sect. 1.2; the objective function f ðxÞ and fitness function F ðxÞ of the new population are then evaluated (rows 11 and 12) according to Algorithms 2 and 1. The twenty-four G-Suite benchmark functions proposed for the IEEE-CEC 2006 competition on real parameter constrained optimization (Liang et al. 2006); A suite of five benchmark mechanical engineering problems (Gandomi and Yang 2011): pressure vessel (Sandgren 1988), welded beam (Deb 1991), cantilever beam (Erbatur et al. 2000), speed reducer and spring design (Gandomi and Yang 2011). The main parameters related to these latter problems are presented in Table 2, where n is the number of decision variables, LI is the number of linear inequality constraints, NI is the number of nonlinear inequality constraints. LE and NE are the number of linear and nonlinear equality constraints, respectively. The complete formulation and mathematical definition of these problems are already well documented in those references; thus, they will not be reproduced here.
To illustrate the advantages of the MCR constrainthandling technique described in Sect. 3.1, the comparisons presented later in Sects. 3 and 4 will highlight the improved performance of the SADE-kNN-MCR and SADE-kNN-MCR 1 variants when applied to the following problems that present constraints with different orders of magnitude (representative of many actual engineering problems as described in Section 3.1): functions G10 and G21 from CEC2006, the pressure vessel and the welded beam engineering problems. Table 3 shows the lower and upper bounds for the constraint functions of those problems: the differences in scale and magnitude are evident. It is interesting to observe that these bounds are analytically determined only for the G10 and G21 functions, and for the Pressure Vessel problem. For the welded beam problem, which presents severely nonlinear constraint functions, orders of magnitude have been numerically estimated.
Following the suggestion of Liang et al. (2006), the algorithms are tested from the results of twenty-five independent executions performed for each benchmark problem. Also following the specifications of Liang et al. (2006), the tolerances for satisfying the equality constraints, and for determining if the optimal solution x* has been found, are, respectively, 0.0001 and f(x*)f(x) \ 0.0001. The maximum number of function evaluations is set to 500,000 for the G-Suite problems. Regarding the structural engineering problems, 10,000 function evaluations are considered to obtain the results presented in Sect. 3.2 to compare the different SADE-kNN algorithms. Table 4 presents the values of the parameters that the proposed method share in common with the other algorithms being compared.

Parameters of the algorithms
The parameters that are specific for the SADE-kNN-MCR and SADE-kNN-MCR 1 variants have been defined from numerical experiments. As described in Sect. 1.2, the strategy that controls the reinitialization of À10 À1 g 6 10 À1 À10 7 g 7 10 4 An enhanced surrogate-assisted differential evolution for constrained optimization problems 6399 the population to avoid stagnation relies on a user-defined parameter that is the number of generations in which stagnation is observed before triggering the population ''rebirth''. The results of the experiments indicated that an appropriate value for this parameter would be 1,000 generations.
Regarding the merit-based database management procedure described in Sect. 1.1, similar experiments have also indicated that an appropriate value for the percentage of the population (%N) comprising the lower-score solutions to be replaced by the exact evaluation of their objective function and introduction in the database at each generation is 22%.
It is interesting to mention that the MCR constrainthandling technique described in Sect. 3.1 has no adjustable parameters.

Performance assessment
The performance of the methods is assessed by direct comparisons and by nonparametric statistical tests, as described next.

Direct comparisons
Direct comparisons are performed by observing the most relevant performance measures collected in tables: best, average, worst results for 25 runs and average of the number of function evaluations to find a feasible solution, considering only the successful runs (succ_nfe); ratio between the number runs in which the algorithms find the optimum solution and the total number of runs (succRate); average of the number of function evaluations (nfe) (when succRate = 1, then succ_nfe = nfe); ratio of runs where the algorithm finds a feasible solution (feasRate).

Sign test
The sign test (ST) is a statistical procedure that indicates if the results of two methods are statistically different. It employs pairwise comparisons that count the number of cases on which a method wins. The method that provides a better result for a given problem receives a positive sign, while the other receives a negative sign. The significance of the difference between these pairwise comparisons (or significance level) is tested by calculating probability values (p-values) (Gibbons and Chakraborti 2003) from the totalized signs. If the p-value is higher than a given threshold value a, then the null hypothesis (H 0 ) that both methods are equivalent is accepted; otherwise H 0 is rejected, and the performance of the methods is different. The most commonly used threshold values are 0.01 and 0.05 (Fisher and Theory of statistical estimation. 1925;Corder and Foreman 2009;Derrac et al. 2011;Lynn and Suganthan 2015), representing a significance level of, respectively, 1% and 5%, i.e., the chance of rejecting the null hypothesis when in fact it is correct.

Performance profiles
The comparisons also use another nonparametric statistical test, specifically the performance profiles (PPs) that were introduced in Dolan and Moré (2002) to evaluate and compare the performance of a set of S solvers for a given set of optimization problems P. The PPs are a useful tool for the easy visualization and interpretation of results from computational experiments. Briefly speaking, they comprise graphs with several curves (one for each solver s [ S considering all problems p [ P), built considering a given metric m p,s . This metric can be the number of evaluations of the objective function; or a statistical parameter such as mean or best value computed from objective function values obtained in several independent runs of the method. The performance ratio r p , s, is then calculated dividing this metric m p,s by the corresponding value of the metric obtained by the solver that performed best on problem p. An overall assessment of the performance of solver s may then be given by: where n p is the total number of problems in set P and q s ðsÞ [ < ? [0, 1] may be seen as the cumulative distribution function for the performance ratio r p , s ; that is, the probability of solver s to have its ratio r p , s within a factor s of the best possible ratio. For s = 1, q s ð1Þ indicates the accuracy of solver s, being a normalized measure of the number of times it was able to provide the best value for the considered metric, or the probability that it will win over all other solvers in set S; q s ð1Þ = 1 indicates that this solver wins over all other solvers in all problems P.
For larger values of s, q s ðsÞ indicates the robustness of the solver, i.e., the number of times it is able to provide feasible solutions; considering an upper limit r M as an arbitrarily large value assigned to r p , s when solver s is not able to find a solution for problem p, then q s ðr M Þ = 1 and the probability that method s solves a problem is given by Overall, the highest the q s values, the better the method is, so a quick visual assessment of the performance of all methods may be performed simply by observing which one has the curves above all others. Conversely, in some cases a more precise numerical indication may be provided by calculating the areas under the curves and then drawing a An enhanced surrogate-assisted differential evolution for constrained optimization problems 6401   It can be seen in Table 5 that the SADE-kNN-MCR 1 method was able to find the optimal solution for 21 out of the 22 functions. Although it had not reached the known optimum for problem G19, the error is relatively small. Also, the fact that the best and worst solutions are very close (that is, the results found in the twenty-five runs present a small variation) indicate the robustness of the method. Considering these 21 functions for which the SADE-kNN-MCR 1 was able to reach the optimal solution, it frequently (except for function G17) obtained the highest success rate (succRate), indicating its best performance in terms of the greater number of runs in which it can find the optimum.
The other methods were able to reach the global optimum for only 15 out of these 21 functions (failing in functions G2, G3, G7, G9, G10, G19 and G21). Considering those 15 functions for which all methods found the optimal solution, the highest succRate (always equal to 1) is again obtained by the SADE-kNN-MCR 1 . The feasibility rate (feasRate) obtained by the SADE-kNN-MCR 1 is again always equal to 1.
Considering the 6 problems for which only SADE-kNN-MCR 1 was able to find the optimal solution, in three of them (G7, G9 and G10), it obtained a success rate of 1 (that is, all runs reached the optimal solution).
Regarding the G10 and G21 problems that present constraints with different units and/or different orders of magnitude, shown in Table 3: for function G10, the enhanced performance of the SADE-kNN-MCR method can be measured by the better average and worst metrics compared with the original SADE-kNN; then, the association of the ''rebirth'' strategy makes the SADE-kNN-MCR 1 obtain the optimal solution in all runs. For function G21, the enhanced performance of the SADE-kNN- MCR and SADE-kNN-MCR 1 methods is reflected by their higher number of feasible runs (FeasRate). Besides the increased performance in finding optimal solutions as commented above, the SADE-kNN-MCR 1 also reduces the average number of function evaluations to find a feasible solution (succ_nfe). To assess this behavior, in the cells corresponding to succ_nfe for SADE-kNN-MCR 1 and considering only the 15 test problems for which all three methods reached at least one successful run, Table 5 shows percentage values in parenthesis that correspond to the ratio between the succ_nfe values provided by SADE-kNN-MCR 1 , and the smaller succ_nfe value observed for the other methods. The number of evaluations is reduced for 13 out of those 15 problems. For the other two problems (G13 and G17), although requiring a higher number of evaluations, SADE-kNN-MCR 1 presents better success rates; for the G13 function, for instance, it finds the optimal solution in 100% of the runs, while SADE-kNN succeeds in only 28%. Table 6 summarizes the SuccRate values for all methods, considering the 15 test problems for which all three methods reached at least one successful run. For all those problems (except G17) SADE-kNN-MCR 1 found the optimal solution for 100% of the runs, while the other methods reached SuccRate = 1 for only 6 out of those 15 problems. The results of the Sign Test from the pairwise comparison between SADE-kNN-MCR 1 and each other method are shown in the last three lines of Table 6 (discarding the draws). The number of wins and losses confirms the observations from Table 5, i.e., the SADE-kNN-MCR 1 outperforms all other methods. All p-values (computed according to the procedure described in Sect. 2.3.2) are small enough to reject the null hypothesis, meaning the pairwise compared samples are not statistically equivalent in all cases. This demonstrates that SADE-kNN-MCR 1 outperforms all other methods.

Performance profiles
The SuccRate values presented in Table 6 were also used as the metric for the generation of the performance profiles, according to the procedure described in Sect. 2.3.3. The PPs shown in Fig. 1 again indicate that SADE-kNN-MCR 1 is the best performer in terms of accuracy and

Fig. 1 Performance profiles-SADE-kNN methods
An enhanced surrogate-assisted differential evolution for constrained optimization problems 6405 robustness: it is easily seen that its curve is above all others, reaching q s ð Þ ¼ 1 with the lowest s value. Regarding the performance of SADE-kNN-MCR compared with the original SADE-kNN method: the contribution of the MCR procedure has already been indicated at some of the direct comparisons of Table 5 (see, for instance, the remarks regarding the G10 and G21 problems in the second-to-last paragraph of Sect. 3.1.1). Also, Table 6 indicates that SADE-kNN-MCR wins for 6 of the benchmark problems (G5, G11, G13, G15, G18 and G23) and loses for 3 problems (G14, G16 and G17). Finally, this better performance of the SADE-kNN-MCR method can also be visualized when comparing the respective PP curves of Fig. 1: the curve corresponding to SADE-kNN-MCR is above the SADE-kNN curve for most of the range of s values; also, it reaches q s ð Þ ¼ 1 with the lowest s value. Table 7 compares the results of the methods for the suite of five benchmark mechanical engineering problems shown in Table 2. All algorithms provided feasible solutions in all runs. Values in boldface indicate the best results among the algorithms.

Suite of structural engineering problems
Regarding the ''best'' solutions, both SADE-kNN-MCR 1 and SADE-kNN-MCR are the best performers, except for the speed reducer problem where all of them achieved the same solution.
It is important to stress that, for the pressure vessel and welded beam problems that present constraints with different units and/or different orders of magnitude, the better performance of the MCR-based methods compared with the SADE-kNN algorithm may be due to the increased efficiency of the MCR constraint-handling method in dealing with this type of problems, as indicated in Garcia et al. 2017b. The increased robustness of the SADE-kNN-MCR 1 compared with the other methods is demonstrated by the average and worst results, since it presents the best performance for all problems except the tension/compression spring.
Here, since the number of examples (only five benchmark problems) is low, statistical performance assessment methods such as the sign test cannot be adequately employed. Even so, it may be illustrative to compare the performance profiles for the methods. Since SADE-kNN-MCR 1 and SADE-kNN-MCR found the same best values for all problems, the average was employed as the comparison metric to be evaluated by the performance profiles. Figure 2 shows that SADE-kNN-MCR 1 finds best values for 80% of the problems with q 1 ð Þ ¼ 0:8.  An enhanced surrogate-assisted differential evolution for constrained optimization problems 6407 Altogether, SADE-kNN-MCR 1 is the best overall performer.  Table 8 compares the results of the methods, in terms of the available main statistical parameters that characterize their performance: the best solution found among the runs; the average number of function evaluations considering the successful runs succ_nfe; and the ratio between successful and total runs SuccRate. The cells corresponding to succ_nfe show also, in parenthesis, the percentage values corresponding to the reduction/increase of the succ_nfe values provided by SADE-kNN-MCR 1 in relation to the smaller value observed for the other methods.

Comparison with other recent methods
In Table 8, it is interesting to compare the performance between the SA-DECV and the SA-DECV-SR methods. The latter is a variation that uses a rank-based constrainthandling technique, instead of simple feasibility rules; thus, unlike the SA-DECV, the SA-DECV-SR provides optimal solutions to problems G7, G10 and G19. It also improves the solution found by SA-DECV in problem G21. SA-DECV-SR also produced a slight improvement in the number of optimal executions. However, the performance of these algorithms is still inferior to those presented by SADE-kNN-MCR 1 and e-DE: SA-DECV methods failed in functions G13, G14, G17 and G23, while SADE-kNN-MCR 1 and e-DE found optimal solution for most problems. Moreover, they found good solutions in most executions, as indicated by their SuccRate. Considering the problems that present constraints with different orders of magnitude, such as G10 and G21 problems, SADE-kNN-MCR 1 presents an enhanced performance. For function G10, besides finding the optimal solution in all runs, it also requires the smaller number of function evaluations. Table 9 summarizes the main results of Table 8, collecting the number of functions for which each method was able to find the optimal solutions in at least one run, and required the lower number of function evaluations to find the optimal solution.
A more comprehensive comparison of the performance of all methods will be presented in terms of nonparametric statistical tests previously described in Sect. 2.3: the sign test and the performance profiles. Table 10 presents the results of the sign test based on the pairwise comparison between SADE-kNN-MCR 1 and each other method in terms of the average number of function evaluations (succ_nfe). The number of wins and losses confirms the observations based on Table 8, i.e., the SADE-kNN-MCR 1 outperforms all other methods. The statistical assessment in terms of the p-values indicates that all of them are below the threshold, which confirms that the superiority of the SADE-kNN-MCR 1 is statistically significant.

Sign test
Regarding the success rate (succRate), both SADE-kNN-MCR 1 and e-DE reach 100% of optimal executions in 16 problems. e-DE gets more optimal runs in 4 out of the remaining 6. However, this leads to a p-value of 0.3437, not small enough to indicate that there is a statistically significant difference between these two methods regarding this particular metric.

Performance profiles
The succ_nfe values were also used as the metric for the generation of the performance profiles shown in Fig. 3, according to the procedure described in Sect. 2.3.3. Differently from the PP curves compared in Sects. 3.1.3 and 3.2. for some curves, it is not quite obvious if they are above or below the other competitors. Therefore, to provide a better global indication of the relative performance of the methods, in this case it is useful to build the bar chart of Fig. 4 indicating the areas under the PP curves.
In Fig. 3 it is interesting to observe that, for smaller values of s, the curve corresponding to the e-DE method is located below those for the SA-DECV-SR and ASRES algorithms; however, since e-DE is able to find feasible solutions for more problems than SA-DECV, SA-DECV-  An enhanced surrogate-assisted differential evolution for constrained optimization problems 6409 SR and ASRES, it ends up with the second larger area in the bar chart of Fig. 4. The overall best performer is SADE-kNN-MCR 1 , since its curve approaches q s ð Þ ¼ 1 with the lowest s value and presents the larger area in the bar chart of Fig. 4, reflecting the fact that it requires fewer evaluations for the greater number of problems.

Suite of structural engineering problems
Finally, this section compares the results for some of the benchmark engineering problems shown in Table 2. Here the proposed SADE-kNN-MCR 1 method is compared only with SA-DECV-SR, since results for the other methods could not be found in the literature. Anyway, this is a useful comparison since both methods share the following characteristics: (1) DE is used as the search algorithm (although with different recombination variants: SADE-kNN-MCR 1 uses DE/rand/1/bin, while SA-DECV-SR combines other DE variants); (2) some of the solutions are approximated by a similarity-based surrogate model (both methods use k-NN, but SA-DECV-SR uses the surrogate model also to approximate the sum of constraint violations); and 3) rank-based constraint-handling techniques are employed (SADE-kNN-MCR 1 uses MCR, while SA-DECV-SR employs the Stochastic Ranking (SR) method (Runarsson and Yao 2000)).
For these comparisons, the maximum number of function evaluations is set to 30,000 function evaluations for the welded beam problem; 20,000 for the pressure vessel; and 5,000 for the tension/compression spring problem. Table 11 compares the statistical analysis of the results for the three benchmark problems considering 25 runs: the best, average and worst solutions among all runs; and the average number of function evaluations considering the successful runs. Unlike Sect. 5.2, where optimal reference values were not used as a stopping criterion, the results presented in parentheses were considered as optimal in order to compare the average number of function evaluations that each algorithm required to reach it.
For the welded beam problem, SADE-kNN-MCR 1 provided the better solution as indicated in the ''best'' column of Table 11. It also obtained this optimal solution in all runs, as indicated by the worst value, which is equal to the best one. Also, its average number of function evaluations (3,983.0) is nearly half of the value required by the SA-DECV-SR.
Considering the pressure vessel problem, SADE-kNN-MCR 1 also provided better solutions regarding all statistical measures. It obtained the accuracy required for the optimal solution in 20 out of the 25 runs, employing an average number of 6,719.7 function evaluations. In (Miranda-Varela and Mezura-Montes 2018) it is not clear whether the value of 6,139 evaluations refers to the number of function evaluations considering only the optimal runs; anyway, it can be concluded that SADE-kNN-MCR 1 is the most efficient algorithm in obtaining the best solutions considering all runs.
Finally, for the tension/compression spring problem, SADE-kNN-MCR 1 provided better results regarding most of the measures shown in Table 11: average, worst, and number of function evaluations. SA-DECV-SR performed better only for best solution, where the value provided by SADE-kNN-MCR 1 (0.012732) reached the first termination criterion (f(x*) -f(x) \ 0.0001), thus stopping the evolution process.
Again, the number of problems is low; thus, the sign test cannot be adequately employed. However, it may be illustrative to compare the performance profiles for the SADE-kNN-MCR 1 and SA-DECV_SR algorithms. Figure 5 shows the performance using as metric the succ_nfe values presented in Table 11, indicating that SADE-kNN-MCR 1 is the best performer, since not only it gives better results for two problems (q 1 ð Þ ¼ 0:666), but still finds optimal executions in all of them (there is s such that q s ð Þ ¼ 1).

Conclusions and extensions
This work presented a new merit-based surrogate-assisted evolutionary method designed for applications to complex real-world constrained optimization problems. Besides presenting a very good accuracy in finding optimal solutions, the main goal of the method is to reduce the number of objective function evaluations. Those goals of accuracy and computational cost savings are pursued by the following approaches: (a) employing the multiple constraint ranking (MCR) technique to deal with the large number of constraints with different orders of magnitude and/or different units that are common in the aforementioned complex engineering problems; (b) devising an enhanced approximation model based on the knearest-neighbor (k-NN) algorithm associated to a meritbased database management procedure; and c) implementing a ''rebirth'' strategy that triggers a random reinitialization of the population whenever it stagnates on local optima.
These approaches have the potential to present significant cost savings by reducing the number of objective function evaluations: the MCR does not require evaluations when there is no feasible solution in the population; the merit-based database management scheme of the k-NN model evaluates only part of the population; and even the population ''rebirth'' strategy does not require additional computational costs, since the new candidate solutions generated whenever the population converges to a local optimum are used to approximate their objective function values.
Two variants of the proposed method (SADE-kNN-MCR and SADE-kNN-MCR 1) were applied to a set of 24 constrained benchmark problems and 5 mechanical engineering problems. The results from both variants were compared with those from the original SADE-kNN method (Garcia et al. 2017a). Then, the best performer among the SADE methods (SADE-kNN-MCR 1) was compared with other algorithms presented in the literature that also employ surrogate models and constraint-handling techniques (SADE-kNN, ASRES, SA-DECV, SA-DECV_SR and e-DE). The performance of all methods is assessed based on direct comparisons, nonparametric statistical tests, and performance profile (PP) graphs. The direct comparisons consider statistical parameters collected from all runs, taking different measures (best, average, worst results, number of evaluations, success ratio, etc.). The nonparametric statistical tests consider the well-known sign test that indicates if the results of two methods are statistically different. Finally, the PP graphs allow a quick global evaluation of the accuracy and robustness of the methods. The comparisons have shown that both variants of the proposed method presented superior accuracy and drastically reduced the number of function evaluations. The best performer is SADE-kNN-MCR 1 since in several executions it was able to escape from local optima.
This method may be particularly suited for real-life engineering problems such as, for instance, the design of offshore systems that have been addressed in Vieira et al. (2012);Lucena et al. 2014;Monteiro et al. 2016;Lima et al. 2005;Pina et al. 2011). Those are complex engineering problems defined on high-dimensional spaces, with many design variables, nonlinear objective and constraint functions (the latter comprised by constraints with different orders of magnitude and/or different units), and requiring evaluations by expensive Finite Element nonlinear dynamic analyses; thus, reducing the number of evaluations is a critical issue. The potential of the SADE-kNN-MCR 1 method to solve such practical engineering problems is indicated by its performance in the welded beam and pressure vessel problems, which are also characterized by severely nonlinear constraint functions with diverse magnitudes. Thus, extensions of this work will incorporate this method into optimal design procedures for mooring systems, submarine pipelines and risers that have been considered in Vieira et al. (2012); Monteiro et al. 2016; Baioco and de Lima Jr. MHA, Albrecht CH, de Lima BSLP, Jacob BP, Rocha DM2018). It is expected that the resulting procedures will lead to increased efficiency in terms of accuracy, computational costs, and more importantly to better engineering design solutions.  An enhanced surrogate-assisted differential evolution for constrained optimization problems 6411 Data availability Enquiries about data availability should be directed to the authors.

Declarations
Conflict of interest R de PG declares that he has no conflict of interest. Beatriz Souza Leite Pires de Lima declares that she has no conflict of interest. ACde CL declares that he has no conflict of interest. Breno Pinheiro Jacob declares that he has no conflict of interest.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.