Voronoi-based Efficient Surrogate-assisted Evolutionary Algorithm for Very Expensive Problems

Very expensive problems are very common in practical system that one fitness evaluation costs several hours or even days. Surrogate assisted evolutionary algorithms (SAEAs) have been widely used to solve this crucial problem in the past decades. However, most studied SAEAs focus on solving problems with a budget of at least ten times of the dimension of problems which is unacceptable in many very expensive real-world problems. In this paper, we employ Voronoi diagram to boost the performance of SAEAs and propose a novel framework named Voronoi-based efficient surrogate assisted evolutionary algorithm (VESAEA) for very expensive problems, in which the optimization budget, in terms of fitness evaluations, is only 5 times of the problem's dimension. In the proposed framework, the Voronoi diagram divides the whole search space into several subspace and then the local search is operated in some potentially better subspace. Additionally, in order to trade off the exploration and exploitation, the framework involves a global search stage developed by combining leave-one-out cross-validation and radial basis function surrogate model. A performance selector is designed to switch the search dynamically and automatically between the global and local search stages. The empirical results on a variety of benchmark problems demonstrate that the proposed framework significantly outperforms several state-of-art algorithms with extremely limited fitness evaluations. Besides, the efficacy of Voronoi-diagram is furtherly analyzed, and the results show its potential to optimize very expensive problems.


INTRODUCTION
The evaluation of solutions to some real-world problems could be difficult or very expensive, in terms of computational cost or money. For instance, evaluating once a solution to the Navier-Stokes equations involves several hours of computational fluid dynamic simulation [1]. Therefore, it is often the case that only a limited amount of budget, in terms of the number of evaluations, is available for searching for a solution. How to search efficiently toward an optimal solution within a small number of evaluations is a crucial problem.
In the past decades, evolutionary algorithms (EAs) have been widely applied to many real-world problems [2]. However, EAs' ability of converging close to the global optimum relies on the availability of a sufficiently large number of Hao Tong, Changwu Huang, Jialin Liu  fitness function evaluations. Obviously, the common EAs is not suitable to solve very expensive problems due to the large number of fitness evaluations (FEs) required to obtain an acceptable solution. Hence, surrogate-assisted evolutionary algorithms (SAEAs), which sometimes sample form a cheap surrogate model instead of successively calling the actual complex evaluation process in conventional EAs, have been proposed to solve computationally expensive problems (CEPs) in mid-1980s [3].
SAEAs have been popular over the past years for its lower number of FEs required for convergence thanks to the use of surrogate models. Many effective models have been employed to reduce the number of FEs required. Examples include polynomial regression (PR) [4], Kriging model [5], support vector machine (SVM) [6], radial basis function (RBF) [7] and neural network (NN) [8]. Besides, the model management strategies play a significantly important role in guaranteeing the convergence of SAEA [9]. One of the most popular strategies is the individual-based model management [10], in which only a few individuals are evaluated by actual function in each generation. Generally, only the individual considered to be the best by the surrogate will be evaluated to exploit the search space [8], but it is likely to be trapped in a local optimum and easily cause premature convergence. On the other hand,information of uncertainty provided by Kriging model is employed to construct much better re-evaluation strategy, like Expected Improvements(EI), Lower Confidence Bound(LCB), Possibility of Improvements(PoI) [11], which are improve the accuracy of surrogate as well as the ability of exploring the search space.
Although SAEASs require fewer FEs to converge than common EAs and variations of SAESs have shown their strength in solving some real-world CEPs (e.g., engineering design [12], health services [13], interactive design [14], etc.), the FEs required is still problematic in some very expensive problems. Take the aircraft design as an example, the duration of accomplishing one crash simulation varies from 36 to 160 hours [15]. If the SAES starts to converge to a global optimal solution after 100 FEs, the total optimization time will be between 150 and 660 days. Despite the importance and demand of handling the budget issue, there is a lack of work considering very limited budget in the literature. For example, Lu et al. [16] applied differential evolution (DE) assisted by rank-SVM for expensive problems given 100D FEs as budget, where D denotes the problem dimension.
The Gaussian process assisted EA proposed by Liu et al. [5] consumed fewer, but still, 50D FEs. More recently, the particle swarm optimization (PSO) assisted by semi-supervised learning (SSL-PSO) [17] and the active learning based PSO [18] were both given only 11D FEs as optimization budget.
All of the existed works use a budget of at least 11D FEs probably due to the minimal number of FEs required to initialize the corresponding algorithm, otherwise these algorithms are not applicable any more. In this work, we consider the extremely expensive problems with only severely limited actual fitness evaluations available and propose a new optimization algorithm named Voronoi-based efficient surrogate-assisted evolutionary algorithm. This work provides the following novel contributions: • To the best of our knowledge, this work is the first to solve very expensive black-box optimization problem with only 5D FEs. With such a limited computational budget, the preference region and model accuracy of corresponding area are much important that we should devote the limited resource to potential optimal area. • We proposed an efficient algorithm for very expensive optimization problem. Leave-one-out cross validation is employed to detect the uncertain area where the current surrogate model is not able to describe much precisely. Meanwhile, the Voronoi diagram is used to partition the search space and force the algorithm to exploit a relatively better area, which makes the algorithm more efficient for very expensive problems with limited fitness evaluations. • The numerical experiment for the proposed framework manifests its superior performance in very expensive black-box problem compared with the state-of-art algorithms. And we also analysed the influence of Voronoi diagram to the framework through empirical results. The remainder of paper is structured as follows. The literature review for SAEAs and the motivation for this work is introduced in Section II. Section III describes the proposed framework and gives some further discussion on the new algorithm. Experiments are presented in IV for comparing the proposed algorithm with the state-of-art algorithms on a set of benchmark problems. Finally, a brief conclusion and the discussion of future work are drawn in Section V.

LITERATURE REVIEW
The conventional evolutionary algorithm is not ideal to tackle CEPs due to the huge fitness evaluation consumption. SAEAs, replacing the actual expensive fitness evaluation by cheap surrogate evaluation, are very effective for CEPS [19] which could categorized into two kinds of SAEAs in the literature, reducing the computational cost in different levels.
a) Category A: The popular framework of SAEAs is inherited from the canonical evolutionary algorithm, which also contains mutation, crossover and selection operator [9]. The surrogate model applied in the framework mainly plays two roles, the approximation and the preselection.
For the fitness approximation, surrogate models mainly influence the selection process that most individuals are evalu-ated by metamodels while remainder individuals are evaluated by the actual fitness function to guarantee the convergence of optimization. For instance, the classification-and regressionassisted differential evolution (CARDE) makes use of classification and regression methods to assist the selection process in differential evolution algorithm [20]. Besides, some algorithms may cluster the population into several clusters and only evaluated the representative individual of each cluster and the other individuals are approximated by the surrogate model [21].
On the other hand, SAEAs employ the surrogate model for preselection, which reduces the computational cost by improving the quality of individuals of the selection process [22]. During the mutation and crossover processes, the framework generates plenty of offspring and the surrogate model is used to preselect the potentially better individuals for selection operator. For example, the differential evolution algorithm employed the preselection strategy (GPEME) to solve the medium scale expensive problems, where only the most promising solution will be re-evaluated by the actual fitness function in every generation [5].
Although the above framework is very effective for expensive problems, the fitness evaluation consumption is still very huge for very expensive problems because the involved population-based algorithm has to evaluate a certain number of individuals to guarantee the convergence and diversity of the algorithm. For example, the CARDE costs 10000 FEs for 30 dimension problems and the GPEME still uses 1000 FEs for problems of 20 and 30 variables. The number of fitness evaluations required to initialize for this kind of SAEA might not be afforded in a very expensive problem. b) Category B: Different from the framework discussed above, another SAEA framework is inherited from the efficient global optimization (EGO) [23]. This category of framework focuses on sampling solutions in the whole search space to improve the accuracy of surrogate model or the quality of the re-evaluated solution, which aims to balance the exploitation and exploration in the optimization process. For example, the classical EGO algorithm uses EI criterion to re-evaluate solutions, which considers the uncertainty of the current surrogate model and the quality of the re-evaluated solution simultaneously.
Recently, researchers proposed different criteria to describe the model's uncertainty. For example, Wang et al. [18] employed the query by committee active learning to assist this kind of SAEA (CALSAPSO), where the solution with the biggest disagreement among various models is regarded as the most uncertain solution. And the local search focuses on the surrogate model constructed by the top 10% evaluated solutions to exploit the current best area.
Although the optimization ability of this category of SAEAs is significantly inferior to the first category of SAEAs, the fitness evaluation is dramatically reduced to obtained an acceptable solution. For example, the CALSAPSO only consumed 11D FEs to solve the expensive problem. Obviously, this kind of SAEAs is more appropriate for the expensive problem.
In this paper, we focus on very expensive problems that the number of fitness evaluations is only setting as 5D in this paper. Obviously, the SAEA based on EGO framework is more appropriate for this problem. However, CALSAPSO, as the most efficient algorithm in the literature is hard to handle CEPs with extremely limited fitness evaluations. The ensemble model to determine the most uncertain solution could only help to improve the accuracy of the surrogate model. As a sequence, the number of fitness evaluation for exploitation will be very small that the performance for very expensive problem will be hindered. Furtherly, the CALSAPSO collects the top 10% best evaluated solutions to build a surrogate model for local search, which only focuses on the region that the selected solutions determine. It has a big problem that the global optimum is likely to locate outside the local area and then the local search is not able to find a better solution for this case.
To address the limitation of CALSAPSO, we employ the leave-one-out cross-validation (LOOCV) and Voronoi diagram to embed in CALSAPSO's backbone, in which the LOOCV can detect the uncertain area of landscape and obtain a better solution at the same time and the Voronoi diagram can guarantee a better solution to be covered by some best Voronoi cells. As a sequence, the new framework is more efficient for very expensive problems. The details of the proposed framework will be introduced in the next section.

Overall framework
According to the previous discussion, we propose a novel model management strategy to assist evolutionary algorithms, named as Voronoi-based efficient surrogate-assisted evolutionary algorithm (VESAEA). The overall procedure of VESAEA is presented in Fig. 1.
We embed LOOCV error surrogate model and Voronoi diagram into the framework, which are much more efficient than the original framework for very expensive problems. The proposed algorithm also contains the two search stages as the global search and the local search. The global search is made up of two stages: one employs the LOOCV error to reevaluate the most valuable individual and another optimizes the surrogate model to select the potentially best individual. On the other hand, the local search applies Voronoi diagram assisted strategy to select the best sample in the local area. Two search processes are controlled by a performance selector which switches the search stage to another one if the current search has no improvement in the previous generation.
In the framework, particle swarm optimization [24] is employed as the optimization algorithm in the global search process. The surrogate used for both global and local search stages in the framework is the RBF model [25].
The proposed algorithm initializes a dataset by Latin hypercube sampling (LHS) [26], at first. Due to the limited fitness evaluations, we use part of fitness evaluations for initialization, which is set as 2D samples in this work. And then the algorithm operates the global search with evaluating two samples on the basis of LOOCV error and RBF surrogate model. The Voronoi-based local search will be applied if the global search makes no improvement during the optimization process according to the performance selector and vice versa. Finally, the algorithm will terminate if it runs out the total FEs setting as 5D in our experiments.
In the following subsections, the detail of two different stages will be described including global search with LOOCV and Voronoi for local search.

Global search with LOOCV a) Construct model with LOOCV error:
The LOOCV is able to measure the sensitivity of landscape and points with high LOOCV error usually located at the rugged area [27]. Therefore, re-evaluating samples with high LOOCV error will not only help to improve the model accuracy but also might discover the local optimal located in the rugged region.
We could easily calculate the LOOCV error for evaluated samples as presented in Eq. (1), whereŷ i is the approximate fitness value according to the surrogate model constructed excluding point (x, y i ).
However, it is impossible to obtain the LOOCV error for unobserved samples in search space. Suggested by [27], we construct a surrogate model by LOOCV error of evaluated samples to approximate the LOOCV error of every point in the design landscape. Assume the dataset: And by leave-one-out experiments, we could obtain the dataset of samples with LOOCV errors as: Sampling with maximal LOOCV error: After completing all leave-one-out cross validation experiments, the surrogate modelê LOO (x) will be built for LOOCV error. The RBF surrogate model with thin plate spline [25] is used both in the cross-validation experiments and LOOCV error model building. Once the model is constructed, the PSO optimizer will be applied to search for the point with the highest approximate LOOCV error in the whole landscape as presented in the Eq. (2) When PSO satisfies the terminal condition, the solution with the highest LOOCV error found by PSO will be re-evaluated by the fitness function and be saved into an archive. c) Construct surrogate model: The second stage in global search is to optimize the global surrogate model followed by individual re-evaluation in LOOCV error model sampling. The global surrogate modelf (x) is still built by the RBF model and optimized by the PSO algorithm. Therefore, the objective function changes into Eq. (3):

d) Select best approximate solution:
Once the PSO completed the optimization process, the best solution found is reevaluated using actual expensive fitness function and then the new sample will be added to the archive.

Algorithm 1: Pseudo Code of Local Search
Input: Samples: S = {xi|i = 1, 2, ..., N } 1 P rand ← Generate a large number of samples in the search space; 2 n ← Dimension of the problem; 3 |P rand | = |S| × n × 1000 4 Each existed point xi constructs one cell Ci; 5 for pr ∈ P rand do 6 pr is assigned to the closest Voronoi cell Ci.
Output: New sample: x * According to the performance selector, the search process will be switched into local search if the global search performs poorly in the current generation. We creatively applied Voronoi diagram approach to promote the efficiency and performance of local search. In mathematic, the Voronoi diagram is used to partition the space into several small regions as described in Section II. In this framework, we regard the evaluated sample of each Voronoi cell as the representative point and the Voronoi cells with better representative points are regarded as better cells. Then, the local search is operated in better cells. The procedure of Voronoi for local search is presented in Algorithm 1.
a) Partition space with Voronoi diagram: The landscape is firstly partitioned by samples in the sample set. Consider a sample set S = {x 1 , x 2 , . . . , x N }, each point constructs a Voronoi cell V i defined in (4): where, d(x, x ′ ) denotes the Euclidean distance between points x and x ′ [28]. Points of each Voronoi cell are closest to the corresponding evaluated sample of this cell. Thus, the boundary between the cells of any pair of adjacent points is the perpendicular bisector of these two points. Due to the irregular shape of Voronoi cell, it is hard to describe the boundaries with one specific equation. The Monte Carlo (MC) simulation is an alternative method to identify the boundaries approximately [29] which is presented in Lines 1-8 of Algorithm 1 where | · | denotes the size of one sample set. A large amount of samples are randomly sampled in the entire search space, then each sample is assigned to a Voronoi cell according to Eq. (4). b) Identify the better Voronoi cell: Then, the top 10% best cells are selected to form a pool of samples C top as presented in Line 9 in Algorithm 1. Assume the number of 10% best cells and total random samples are k, N P respectively. Then, the C top can be described as Eq. (5).
c) Select best solution in local area: Then, the framework will select the best sample in the set of C top according to the approximate modelf (x) that the objective function is  6). The new sample will also be re-evaluated with expensive fitness function after finishing local search and then added to the archive.
Take Fig. 2 as an example, P 1 , P 2 are two best samples so that C 1 , C 2 are two best Voronoi cells and S is the optimal solution. According to the local search, the new samples will be selected from C top = C 1 ∪ C 2 . As shown in Fig. 2, we could easily find the optimal solution is likely to locate in the few top best Voronoi cells.

Discussion
In order to obtain a better solution, it is very important to trade off the exploration and exploitation. In this framework, a performance selector is applied. The global search explores potentially better area in the whole search space and then the local search is applied to exploit the specific area to obtain a local optimum if global search has no improvement in current generation.
The global search firstly starts with sampling with respect to the LOOCV error to increase the model accuracy. After re-evaluating the sample with maximal approximate LOOCV error, the new surrogate model built by all samples in the dataset is more likely to describe the valuable rugged area of landscape where the local optimal might hidden. Then, reevaluating the best individual of new accurate surrogate model is easy to find a better solution. Therefore, the global search has the ability to not only explore the uncertain area, but also has the potential to obtain a higher quality solution.
By contrast, the local search considers exploiting the potential area. In the framework, the space partition makes local search more focus on the specific better area and meanwhile, the top 10% cells are selected to avoid trapping into a local optimum. Monte Carlo simulation is used to describe the Voronoi cells and it only requires to select the new solution among the generated random samples belonging to the selected cells, which avoids searching next sample by complex algorithms.

NUMERICAL EXPERIMENT AND ANALYSIS
In order to analyse the efficacy of VESAEA, we conduct a series of empirical studies in this section. We compared the proposed algorithm with several state-of-art algorithms on classical some benchmark problems at first. Then the empirical studies of Voronoi diagram's effect to the proposed algorithm are presented.

Experimental Setup
To assess the performance of SAEA, we select one efficient classical algorithm and three state-of-the-art algorithms, which have different characteristics on computationally expensive problems. A brief description of four algorithms is presented below.
• EGO-LCB: Lower confidence bound (LCB) criterion based efficient global optimization [23] is a classical algorithm but performs efficiently on expensive problems, especially on low-dimension problems. • GPEME [5]: DE assisted by Kriging model with individual-based evolution strategy. The LCB criterion is also employed as the re-evaluation strategy. • SSLAPSO [17]: PSO is enhanced by semi-supervised learning which is employed to make use of both evaluated and un-evaluated samples. Two RBF surrogate models are built by true evaluated and approximately evaluated samples respectively to determine new individuals' fitness. • CAL-SAPSO [18]: The committee-based active learning is applied to SAEA framework by using ensemble of PR, RBF and Kriging model. Three surrogate models are used simultaneously to determine the most uncertain samples and potentially best solutions. The comparative experiments are performed on 5 widely used benchmark problems of 5 different dimensions listed in Table II. The size of initial samples is set to 2D and all algorithms terminate after 5D real fitness evaluations. The surrogate models used in all algorithms are all imported from SURROGATES tool-box [30], which is based on a RBF toolbox [31].

Comparative Experiment on Benchmark Problems
The results obtained by five algorithms over 25 independent runs are shown in Table I  where the x-axis ranges from 2D to 5D because the first 2D fitness evaluations are used for initialization, which is same for all algorithms in each run for fair comparison. The average rank is calculated according to the best fitness on each test problem. Then, the Friedman statistical test and its postdoc test, Nemenyi test [32], are both carried out with a 0.05 significance level. The p-value by Friedman test equal to 6.193e − 12 lower than 0.05 significantly, revealing significant difference among five algorithms. Then the Nemenyi test result is presented as adjust p-value in Table I where VESAEA is the control method. It also shows the performance of VESAEA on benchmark problems is greatly different with ECG-LCB,  The Sphere problem is a uni-modal problem and the Griewank problem is a little similar to sphere problem except for many small peaks in the whole landscape. The performance of the proposed framework is also very similar in these two functions whose evolution curves are presented in Fig. 4. Compared with the other four algorithms, VESAEA obtains the best results both in convergence behaviour and the quality of final solution for D = 10, 15, 20, 30. By contrast, the EGO-LCB is slightly more efficient than others when the dimension D = 5. As a sequence, we could generally conclude the proposed model management strategy is efficient in the overall uni-modal whether with or without small peaks.
For Ackley function, which has a nearly flat outer region and a large hole at the centre, the VESAEA is slightly better than other four algorithms for D = 15, 20, 30 while the EGO-LCB performs best for D = 5, 10 as shown in the first column of Fig. 5. Also for D = 15, 20, 30, we could easily observe that CALSAPSO performs well at the beginning of search process but is surpassed by VESAEA in the latter stage of the search, which indicates that the new framework is capable to balance better the exploitation and exploration. Different with the above problems, Rosenbrock function is a multi-modal problem with a very narrow valley. VESAEA obtains the best results for problems of D = 5, 10, 15 and the convergence profile is presented in the second column of Fig. 5. From evolution processes, it is obvious that VESAEA outperforms other four methods totally for problems of D = 10, 15. And for the lowest dimension case, VESAEA is also not very efficient in the early search stage but still gets best result finally similar to former cases in Ackley problems. However for D = 20, 30, the CALSAPSO converges slightly better than VESAEA in the latter stage of the search and obtains higher-quality results.
For the most complicated test problem, Rastrigin function, the performance of VESAEA framework is better than EGO-LCB, GPEME and SSLAPSO and almost equal to CAL-SAPSO of D = 10, 15, 20, 30 as shown in Fig. 6. And when dimension equal to 5, the LCB-EGO still obtains the best solution compared with other four algorithms. Actually, the Rastrigin problem has many significant peaks that have strong capacity of trapping the optimization algorithm in one of them. And all five algorithms are hard to jump out of the local optimal in limited fitness evaluations.
From the above observations, we can see that the VESAEA is capable of solving the uni-modal problems even though with little noise in the search space. VESAEA obtains the best result in 16 out of 20 test cases of D = 10, 15, 20, 30 which indicates its capability in optimizing low-dimension of problems. But the classical framework, EGO-LCB, is always the best algorithm in very low-dimension problems, like D = 5. Besides, experimental results also indicate that it is hard for GPEME and SSLAPSO to find a high quality solution in such limited computational resources. However, the performance of VESAEA deteriorates in very complicated problems with many attractive local optimum, like Rastrigin functions. For this kind of problems, the global search hardly jumps out of the local optimums with limited fitness evaluations and on the other hand, the Voronoi based local search only focuses on the better area in the current generation. As a sequence, the algorithm is easily trapped in the local optimal especially for relatively high dimension cases.

The effect of Voronoi-based local search
In this paper, we creatively employed Voronoi diagram for local search in SAEA framework. In order to study the effect of Voronoi based local search furtherly, we analyzed the average contribution ratio of local search and global search stages over all 25 independent runs. The results are presented in Fig. 3 where the red dotted line represents bisector of contribution ratio.
In Fig. 3, there are 19 out of 25 test instances in which the local search contributes more than the global search, which indicates the potential performance of local search. We can find in Griewank and Sphere problems, the Voronoi-based local search accounts for a large proportion (over than 70%) of the whole search process, while in Rastrigin and Rosenbrock problems, both search stages make about the same contributions to the final result. Moreover, we calculate the relation the between dominated search stage and the optimization result as presented in Table III   To further assess the capability of Voronoi-based local search, we compare VESAEA and VESAEA without Voronoibased local search stage (VESAEA-woVLS) on 5 test problems of dimension d = 15 with 25 independent runs. The results are presented in Table IV and the convergence curves are presented in Fig. 7 where the initial sample size is 2D and the total fitness evaluation is 5D. From the Table IV and Fig. 7, it is very clear that VESAEA performs significantly better than VESAEA-woVLS in Sphere, Ackley and Griewank problems, probably because the better Voronoi cells usually cover the area containing global optimum instead of local optimum. And it has a similar performance to VESAEA-woVLS in Rosenbrock problems because the best Voronoi cells might all fall into the narrow valley and VESAEA can not benefit much from the Voronoi-based local search. But in Rastrigin problem, VESAEA-woVLS performs slightly better than VESAEA, probably because VESAEA-woVLS only focuses on the global search and easily jumps out of many attractive local minimums. In addition, we could find test problems in which VESAEA performs better than VESAEA-woVLS are consistent with problems in Fig. 3 in which the local search contributes more than the global search.
Overall, we could conclude that the Voronoi-based local search is able to improve the convergence of optimization, especially for uni-modal problems and multi-modal problems with many small peaks. And the Voronoi-based local search is not obviously outstanding in complicated multi-modal problems with many attractive local optimums.

CONCLUSION
A Voronoi-based efficient surrogate assisted evolutionary algorithm framework is proposed in this paper for very expensive problems with extremely limited computational resources.     The framework mainly has two search stages including leaveone-out cross validation based global search and Voronoi based local search. The global search stage employs LOOCV and RBF model to improve the accuracy of surrogate model and explore the search space. In local search stage, the Voronoi diagram, which divides the whole space into many cells, is applied to assist the exploitation in local area. We designed a performance selector which switches between two search stages by detecting whether current search stage has improvements or not to trade off the exploitation and exploration. The efficacy of proposed framework is examined by comparing with a few state-of-the-art SAEAs on 25 commonly used test problems. The results demonstrate that the proposed framework is powerful at solving uni-modal problems, even though there are several small peaks around the response surface. Moreover, we evaluates the effect of Voronoi-based local search and the results show that Voronoi diagram contributes a lot in finding a better solution on a very limited computational budget.
The ability of Voronoi for boosting optimization in expensive problems has been demonstrated in this work. In the future, the efficacy of Voronoi-based optimization algorithms will be tested on more complicated test functions. On the other hand, more strategies with Voronoi assisted SAEA for very expensive problems will be considered. Since Voronoi diagram divides the whole search space into many cells, it is highly desirable to investigate more on Voronoi diagram for large-scale problems.