A two-level global optimization method based on hybrid metamodel for expensive problems

Abstract

With the wide application of simulation and optimization tools in engineering problems, how to build a metamodel, which can satisfy the high accuracy requirements of real working conditions and realize fast optimization in the entire design space, becomes a hot issue. On the basis of sequential sampling optimization and updating design space optimization, a two-level global optimization method based on hybrid metamodel is developed. In this method, the hybrid metamodel is constructed using random sampling. In the first level, a space reduction strategy is proposed to reduce the design variable space and guide the search to the promising region. In the second level, Adaptive simulated annealing algorithm is integrated with metamodels to search the global optimal value in the promising region. Several global optimization problems and a real industrial design optimization example are utilized to demonstrate the superior performance of the proposed method.

Introduction

At present, optimization tools such as finite element analysis (FEA) and computational fluid dynamics (CFD) have been widely used in product optimization design. However, the optimization accuracy and efficiency are still the bottleneck for its development. It is reported that it takes Ford Motor Company about 36–160 h to run one crash simulation; ¹ One run of a whole car impact simulation still costs more than 10 h to simulate a 0.1-s impact process. ² It can be seen from the literature that a cheap approximate model is often used as a surrogate of the expensive model to carry out optimization design. ³ Approximate model, which is also known as metamodel, surrogate model or response surface, can relieve the computing budget and improve the design efficiency of simulation-based problem. ⁴ The commonly used approximate models include quadratic response surface model (QRSM),^5,6 kriging (KRG),^7–9 radial basis function (RBF),^8,10,11 and multivariate adaptive regression spline (MARS) ¹² . When the accuracy of metamodel meets the design requirements, the optimization design can be conducted based on the model. In order to establish high accuracy of approximate model with sampling points as less as possible, there are two commonly used methods: the sequential sampling and updating optimization space.

Sequential sampling optimization

Traditional sampling methods can be divided into one-stage sampling and sequential sampling. Because it is difficult for one-stage sampling to obtain suitable sample size, in order to address this issue, researchers seek to improve model accuracy through sequential sampling. ¹³ Currin ¹⁴ carried out sequential sampling based on entropy criterion to maximize the information amount of sample points in the framework of Bayesian method. Osio and Amon ¹⁵ proposed Bayesian surrogate models and a sequential sampling method based on the A-optimality criterion. Lin et al. ¹⁶ proposed the sequential exploratory experimental design (SEED) method, in which data points are identified sequentially based on Bayesian entropy sampling by removing the stationary assumption and introducing correction factors in the calculation of correlations between points. Yao et al. ¹⁷ proposed a gradient-based sequential radial basis function neural network (RBFNN) modeling method, which utilizes the gradient information of the present model to expand the sample set and refine the model sequentially. Kitayama et al. ¹¹ presented a sequential approximate optimization (SAO) procedure, in which the optimum of the response surface and the global minimum of the density function are taken as new sample points to improve the model accuracy. Aute et al. ¹⁸ developed a KRG metamodel for leave-one-out (LOO) error, point with highest estimated value of LOO error is added to the sample set, and maximum entropy design is integrated with maximin distance to improve the model accuracy. Xu et al. ¹³ used Voronoi diagram to partition the design space; the error behavior of each cell is estimated through LOO cross-validation approach; new points are sampled sequentially in a Voronoi cell with the largest error value to rapidly improve the metamodel accuracy.

The sequential sampling methods for single metamodel mentioned above can effectively improve the accuracy of approximate model. However, single metamodel has its specific applicable conditions. In order to exploit the approximate performance, the hybrid metamodel has gradually become a research focus. Zhang et al. ¹⁹ determined the weight of each contributing surrogate model based on the local measure of accuracy for that surrogate model in the pertinent trust region and selected weights for the three component surrogate models adaptively to improve the overall accuracy of the model. Tang et al. ²⁰ constructed hybrid metamodel for black box functions by combining quadratic polynomial response surface with RBF and conducted optimization design by means of particle swarm algorithm. Sun et al. ²¹ considered and compared typical surrogate approaches such as response surface methodology, KRG, RBF, and support vector regression for the best modeling accuracy. Yin et al. ²² integrated the dynamic ensemble metamodeling method with the multi-objective particle swarm optimization (MOPSO) algorithm to improve the accuracy of multi-objective crashworthiness optimization. Fang et al. ²³ established the ensemble of surrogates to approximate the fatigue life function of a truck cab for the purpose of decreasing the modeling error. Ferreira and Serpa ²⁴ created ensemble of metamodels based on least squares (LS) which is used to estimate the ensemble weights without using any explicit error metrics. Gu et al. ²⁵ proposed a hybrid and adaptive metamodel (HAM)-based global optimization method, which applies three metamodels concurrently, sample data points are selected adaptively according to the calculated values of the three metamodels to gradually improve modeling accuracy. Li et al. ²⁶ concentrated on the high-dimensional optimization problems. Multi-surrogate sampling strategy is utilized to construct more reliable model, and the two-level optimization strategy is adopted to improve the efficiency of the surrogate with limited samples. Qi et al. ²⁷ claimed the accuracy of the surrogate model is improved by constructing ensemble surrogates, and the reliability of the solution obtained by the optimization model could be improved using sequential sampling. Ye et al. ²⁸ created an ensemble of surrogates with optimized weight factor to provide robust prediction performance and selected the promising point to supplement the sample set to improve the accuracy of the model.

By adopting sequential sampling method, the accuracy of single metamodel or hybrid metamodel can be improved. For single metamodel, the model accuracy is easily influenced by the characteristics of the approximated object. For example, it is difficult for high-dimensional non-linear function to get high-accuracy model with QRSM, which is more suitable for low-order function. As for sequential sampling method based on hybrid metamodel, difficulty is encountered in trying to obtain proper weight value for each model when setting weight and conducting linear combination. Meanwhile, for sequential sampling method, the clustering situation is easy to appear with the increase in the sample points. Jin et al. ⁵ claimed the accuracy of an approximate model is not decided by the sample size. An interpolation method, like RBF, tends to over-fit the performance, which may lead to large local errors. However, even the large sample size does not necessarily reduce local errors. Therefore, with a certain sample size, reducing and partitioning design variable space becomes another important direction to improve the optimization accuracy and efficiency.

Updating design space optimization

Updating design space optimization is a strategy which gradually updates optimization design variable space and focuses on the key concerned area. Srinivas and Patnaik ²⁹ improved the optimization efficiency of genetic algorithm by reducing variable space. Chen and Smith ³⁰ used representation-level constraints to beneficially reduce the search space and improve the search efficiency of genetic algorithm. Wang et al. ³¹ proposed adaptive response surface method (ARSM), in which the quadratic approximation model is created for computation-intensive design objective function, a threshold or a cutting plane is used to reduce the design space to approach optimal value. In order to address the issue that ARSM requires a great amount of design experiments for high-dimensional design problems, Wang ³² took a general sampling strategy which inherits sample point as long as it falls into the updated design space. Stander and Craig ³³ constructed a series of linear metamodels in a subregion of the design space and iteratively contracted the size of this subregion to approach the optimum. Wang et al. ³⁴ proposed a new global optimization method, in which quadratic regression is performed to detect the key region, and mode-pursuing sampling method is used to generate more sample points in the concerned area. Melo et al. ³⁵ developed domain optimization algorithm (DOA), which can efficiently eliminate search-space regions where there is low probability in containing a global optimum and speed up the search to the global optimum. Wang and Timothy ³⁶ proposed an intuitive methodology which the fuzzy c-means clustering method is utilized to cluster sample points, and optimization design is carried out after systematically reducing the design space to a relatively small region. Zhu et al. ³⁷ applied fuzzy c-means clustering method to identify the reduced attractive regions in the original design space, which accelerates the search of the global optimum. Li et al. ³⁸ proposed a new metamodel-based global optimization method, in which fuzzy c-means method and Gath-Geva clustering method are applied to divide the design space into several small interesting cluster spaces for low- and high-dimensional problems, respectively, the modeling efficiency is improved by means of the reduction principle. Wu et al. ³⁹ divided the entire domain of design variables into multiple cells, and identified rougher cells with large modeling error, and then added a small number of additional sample points into these cells to improve the model accuracy. Long et al. ⁴⁰ integrated ARSM with intelligent space exploration strategy (ISES) to carry out space reduction for higher global optimization efficiency. Gan et al. ² applied two different search strategies, respectively, inside and outside the promising region to improve optimization efficiency. Dong et al. ⁴¹ classified design variable space into global space, reduced medium space, and local space, and the search process is carried out alternately in three spaces to improve search efficiency. Wang et al. ⁴² discussed several novel space reduction methods and their application in sheet metal forming design due to their high efficiency.

By reducing variable space or emphatically searching in the promising region after partitioning the design space, the metamodel-based optimization efficiency can be improved effectively. The optimization method based on space updating can accelerate convergence to the optimal solution. However, the accuracy is influenced by the method of space reduction or space partition.

On the basis of the research on these two methods which aim to improve optimization accuracy and efficiency, this article proposes a hybrid metamodel–based two-level global optimization (HMTGO) method. In the first level, considering the diversity of surrogates, a space reduction strategy based on hybrid metamodel is proposed to reduce variable space and guide the search to the promising region. Compared with the existing space reduction strategy, in which single metamodel and clustering method are utilized, the space reduction strategy proposed in this article can simplify the search process and ensure the reliability of optimization results. In the second level, adaptive simulated annealing (ASA), ⁴³ which can effectively avoid being trapped into the local minima, is integrated with hybrid metamodel to search the global optimal value and ensure the accuracy of optimization results.

Through space reduction strategy and two-level optimization, the method maintains a balance between optimization accuracy and efficiency, and conducts a beneficial research on the engineering optimization problem.

The hybrid metamodel-based two-level global optimization method

In recent years, Latin hypercube design (LHD) sampling is considered as one of the most popular design of computer experiments methods. Optimal Latin hypercube design (OLHD) not only has the projective characteristics of LHD but also has a good space-filling property. ⁴⁴ Therefore, OLHD is adopted to generate sample points in two-level optimization.

Haftka et al. ⁴⁵ claimed that fitting multiple surrogates to the same data is very straightforward in approaches that zoom on promising regions of the design space to restrict the search space. Hence, multiple approximate models are used for space reduction and promising region searching in this article. Two metamodels are not enough to ensure the robust of prediction performance, and more metamodels will increase the amount of computation. According to the test results and the existing literatures, three to five metamodels are suggested. The selected metamodel depends on the type of objective functions and the characteristic of the metamodel itself, and the hybrid metamodel should be suitable for the pre-optimized function. Jin et al. ⁴⁶ pointed out that QRSM is very accurate for low-order non-linear functions and insensitive to the numerical noise; KRG is better for high-order functions; the average accuracy and robustness of RBF are influenced only slightly by sample size. So, these three approximate models can cover quite a few types of functions and are utilized to construct the hybrid metamodel.

ASA is an efficient stochastic global optimization algorithm. This article utilizes this algorithm in the second-level optimization to search the optimum for the following reasons: (1) ASA is not only especially suitable for the optimization of multi-modal, discrete, non-linear, and non-differentiable target functions but also ensures a high search efficiency. (2) Compared with other optimization algorithms, ASA can avoid getting trapped in local optimal effectively. ⁴⁷

Methodological model

The hybrid metamodel-based two-level global optimization method includes two-level search process. As shown in Figure 1, first, conduct data sampling using OLHD, establish three approximate models of QRSM, KRG and RBF, respectively, for time-consuming practical problems. In the first level, evaluate sample points generated from OLHD via these three approximate models; after obtaining the approximate response values, select partial sample points with the lowest response values to reduce variable space. With the continuous reduction of space, the sample points for establishing approximate model are gradually approaching to the true optimal value. In addition, the accuracy of the approximate model is being improved constantly. In the second level, obtain a small variable space after dealing with the reduced space of the first level and conduct optimization search through ASA in the small space by means of three latest approximate models established in the first-level optimization, and then obtain the approximate optimal value of each model. Make a comparison between the approximate optimal value of the latest QRSM figured out by ASA in the small space and the approximate response value of QRSM founded from the space center sample point. If the difference of these two approximate values meets the accuracy requirements, then the searching optimization is terminated; the sample point corresponded to the optimal solution of the latest QRSM figured out by ASA in the small space is the global optimal point; if not, then confirm the updated variable space based on the approximate values founded through ASA on the latest three models. Reduce the space via the latest QRSM model which was constructed in the first-level optimization. Construct updated QRSM in the reduced space and obtain the optimal value via the space reduction strategy again. Repeat the second-level optimization until the optimal value meets the accuracy requirements and then the search process terminates.

Figure 1.

Methodological model of HMTGO.

The specific implementation process

• Step 1. Set the pre-optimized function $f(x)$ and variable space $S$ , generate $Q$ sample points in $S$ using OLHD, and the number of sample points should meet the requirement of constructing the three approximate models at the same time. Meanwhile, in order to carry out LOO cross-validation in step 3, set

$$Q=\frac{(n+1)(n+2)}{2}+1$$

(1)

where n is the number of design variables. Evaluate these sample points using the real function and establish QRSM, KRG, and RBF, respectively.

• Step 2. Continue to generate $M$ sample points using OLHD, and $M$ is obtained according to formula (2)

$$M=\underset{i=1}{\overset{n}{\Pi }}\frac{x_i^u-x_i^l}{k}$$

(2)

where $x_i^u$ and $x_i^l$ are the upper bound and the lower bound of ith design variable. In general, $M$ depends on the dimension of function and the range of design variables. The value of $M$ should not be too large nor too small, the global optimal value is easily missed if the value of $M$ is too small, the computational cost will be increased, and the efficiency will be affected if the value of $M$ is too large. Moreover, according to test experience of a number of functions with different dimensions, the value of $M$ is suggested to be controlled among 50–200. In addition, the boundary difference of each design variable should be divided exactly by $k$ . In order to avoid expensive computational cost of the high-dimensional pre-optimized functions, take $M$ equals 100 for this kind of functions if $k$ cannot meet the above requirements. Then, use the three approximate models to evaluate the sample points, respectively, and obtain the approximate response values.

• Step 3. Rank the approximate response values according to the results received from metamodels in ascending order, then select the first $P=\rho ·M$ response values in QRSM, KRG, and RBF, respectively. $\rho$ can be different values, and the optimization result of different $\rho$ will be compared in section “Function tests.” Confirm the corresponding highest and lowest values of each design variable. Set them as the new upper bound and the lower bound. Calculate generalized mean square cross-validation error (GMSE) of each approximate model and use WTA3 proposed by Goel et al. ⁴⁸ to assign weight to the upper and lower boundary values of different design variables evaluated by the three approximate models. Obtain the new variable space $S\prime$ , and the optimization space $S\prime$ is smaller than $S$ . If the search process moves along at a slow rate thanks to the fact that there is only slight change between $S\prime$ and $S$ , according to the description of optimization space reduction, ³² it should make a 5%–10% reduction on either bound of any of the variables in each iteration to improve the space reduction efficiency. Which means to reduce the selected sample points $P$ from $50\%M$ to $25\%M$ at a rate of 5%, then enter the next iteration.

• Step 4. In $Q$ sample points, inherit the sample points that fall into the new variable space in previous iteration and generate new sample points using OLHD; after evaluating these points with the real function, update approximate models and then repeat step 2 and step 3. When $P<10$ , in order to improve the optimization accuracy, set $M=100$ in step 2 for later iterations and the number of sample points selected in step 3 changes from $50\%M$ to

$$P=\frac{\sqrt[n]{\underset{i=1}{\overset{n}{\Pi }}{x}_i^u-x_i^l}\times 50}{\sqrt[n]{\underset{i=1}{\overset{n}{\Pi }}{x}_{imax}-x_{imin}}}$$

(3)

where $x_{imax}$ and $x_{imin}$ are the initial maximum boundary and the initial minimum boundary of ith design variable. If $P-P\prime \le 5$ , in other words, if the difference of selected sample points in the first two iterations of step 3 is no more than 5 after resetting $M=100$ , then select 10 points when repeating step 3 in later iterations until $P$ is less than 10, subsequently select points according to formula (3).

• Step 5. When $P=0\mathrm{or}1$ , continue the optimization iteration. According to Taylor’s theorem, any smooth function can be locally, accurately approximated by a quadratic function. ³⁴ Meanwhile, Liu et al. ⁴ pointed out the fact that QRSM uses the LS estimation method to ensure the sum of the prediction errors is minimized, so it provides better local accuracy. Because the corresponding variable space at this time is small enough, QRSM always has higher accuracy compared with the other two approximate models. The optimization result of functions in section “Function tests” also verifies the conclusion. So, the obtained optimal value of QRSM is the searching result of the first-level optimization.

• Step 6. Increase and decrease the optimal values of each design variable in step 5 in a certain proportion, then obtain a small variable space. By reference to the literature, ⁴⁹ the variation equals to 2.5% of the entire variable space for each design variable. Specific implementation strategies are as follows:

Set $x_i$ the optimal value for the ith design variable, and the updated space $[x_i^l,x_i^u]$ is

$$x_i^u=x_i+(x_{i\max }-x_{i\min })\times 2.5\%$$

(4)

$$x_i^l=x_i-(x_{i\max }-x_{i\min })\times 2.5\%$$

(5)

In the updated space, on the basis of the approximate models constructed in the latest iteration, take ASA algorithm to search for the optimal value and obtain the optimization result. If it meets the termination criterion, the entire search process terminates; else, go to step 7.

• Step 7. In the previous step, among the optimization results received in the updated space with ASA, with the result of QRSM as a reference, select the optimal value between KRG and RBF which is closer to the result of QRSM and then determine the new variable space based on the optimal values of the selected two approximate models. If any design variable has the same boundary value, then decrease the number of decimal places of boundary value according to the optimization accuracy. Since the number of decimal places is decreased by the premise of ensuring accuracy, it has little impact on optimization results. More importantly, it is convenient for OLHD to generate new sample points. Generate 100 sample points using OLHD in new space, obtain the approximate response values according to the results of QRSM, select 10 sample points with the minimum response values, determine the new upper and lower boundary value of the design variable via these sample points, and obtain the updated variable space.

• Step 8. Construct updated approximate models by sampling data points in the updated variable space in step 7, generate 100 sample points using OLHD. After evaluated by QRSM, the sample point with the minimum approximate response value is chosen as the optimal point for the second-level optimization. Process the optimal value in accordance with step 6, if the optimal value meets the termination criterion, the entire optimization process terminates; else, repeat step 7 and step 8 until the termination criterion is met.

The flowchart of the method is given in Figure 2.

Figure 2.

Procedures of HMTGO.

Termination criterion

Among the optimization results, the approximate response value of QRSM is adopted for reference, when

$$|{\hat{y}}_{Q-\mathrm{ASA}}-{\hat{y}}_Q|\le \epsilon$$

is met, where ${\hat{y}}_{Q-\mathrm{ASA}}$ is the optimal value figured out by ASA algorithm in the updated space $[x_i^l,x_i^u]$ in step 6 based on QRSM model, ${\hat{y}}_Q$ is the approximate response value of QRSM founded from the optimal point. The search process terminates.

Case study

At present, there is not a standard set of test problems for global optimization algorithms, some test functions are widely used.^4,11,25,31 Among those, the Goldstein and Price (GP) function is a high-order polynomial function which is frequently used to test new optimization algorithms for their convergence characteristics. ⁵⁰ This function rises up to higher than ${10}^9$ and has four local minima in a small region. The global optimal point is $x_1=0,x_2=-1$ , and the expression is as follows

$$\begin{array}{l}f(x_1,x_2)=(1+{(x_1+x_2+1)}^2\\ (19-14x_1+3{x_1}^2-14x_2+6x_1{x}_2+3{x_2}^2))\\ ({(30+2x_1-3x_2)}^2\\ (18-32x_1+12{x_1}^2+48x_2-36x_1{x}_2+27{x_2}^2))\\ x_1,x_2\in \left[-2,2\right]\end{array}$$

This article takes one set of test data of GP function as an illustration and explains the space reduction strategy for inheriting sample points and the two-level optimization strategy. First, due to GP function has two design variables; generate seven initial sample points using OLHD according to formula (1) to construct QRSM, KRG, and RBF model; consider the variable space $x_1\in [-2,2],x_2\in [-2,2]$ ; follow formula (2) to generate 64 sample points using OLHD; obtain the approximate response values of these sample points according to the results evaluated with approximate models; rank these values in ascending order; select the first 32 sample points with the lowest response values, respectively; and find the corresponding upper and lower boundary value of design variables. Assign weight to the upper bound and the lower bound of each design variable according to the results received from the three approximate models, obtain the updated variable space (if the updated space cannot satisfy the space reduction requirement, continue updating variable space according to step 3 until the requirement is met, go to next iteration). The specific variable spaces, number of sampled points, and number of reserved points for each iteration are listed in Tables 1 and 2.

Table 1.

GP function optimization data for the first level.

Iteration	Variable space	Number of sampled sample points	Number of reserved sample points
0	$x_1\in [-2,2],x_2\in [-2,2]$	64	32
1	$\begin{array}{c}{x}_1\in [-1.4644701.988190]\hfill \\ x_2\in [-1.9190401.016229]\hfill \end{array}$	41	20
2	$\begin{array}{c}{x}_1\in [-1.4224101.815045]\hfill \\ x_2\in [-1.9190400.462911]\hfill \end{array}$	31	15
3	$\begin{array}{c}{x}_1\in [-1.4224101.815045]\hfill \\ x_2\in [-1.0447600.462911]\hfill \end{array}$	20	10
4	$\begin{array}{c}{x}_1\in [-1.3190701.591440]\hfill \\ x_2\in [-0.8256300.320058]\hfill \end{array}$	100	23
5	$\begin{array}{c}{x}_1\in [-0.6440600.409748]\hfill \\ x_2\in [-0.8256300.207950]\hfill \end{array}$	100	13
6	$\begin{array}{c}{x}_1\in [-0.6412600.293918]\hfill \\ x_2\in [-0.820410-0.030540]\hfill \end{array}$	100	11
7	$\begin{array}{c}{x}_1\in [-0.4071200.094956]\hfill \\ x_2\in [-0.810140-0.668500]\hfill \end{array}$	100	3
8	$\begin{array}{c}{x}_1\in [-0.194560-0.113900]\hfill \\ x_2\in [-0.807980-0.796150]\hfill \end{array}$	100	1
9	$x_1=-0.144050,x_2=-0.807620$	N/A	N/A

GP: Goldstein and Price.

Table 2.

GP function optimization data for the second level.

Iteration	Initial point	Processed small variable space	Updated variable space	Optimal point
10	$[-0.14405,-0.80762]$	$\begin{array}{c}{x}_1\in [-0.244050.04405]\hfill \\ x_2\in [-0.9076-0.70762]\hfill \end{array}$	$\begin{array}{c}{x}_1\in [-0.1461-0.04405]\hfill \\ x_2\in [-0.9076-0.81787]\hfill \end{array}$	$[-0.05278,-0.90653]$
11	$[-0.05278,-0.90653]$	$\begin{array}{c}{x}_1\in [-0.152780.04722]\hfill \\ x_2\in [-1.0065-0.80653]\hfill \end{array}$	$\begin{array}{c}{x}_1\in [-0.028180.04722]\hfill \\ x_2\in [-1.00315-0.9306]\hfill \end{array}$	$[0.02384,-0.99227]$
12	$[0.02384,-0.99227]$	$\begin{array}{c}{x}_1\in [-0.076160.12384]\hfill \\ x_2\in [-1.0922-0.89227]\hfill \end{array}$	$\begin{array}{c}{x}_1\in [0.0014730.024004]\hfill \\ x_2\in [-1.0009-0.97265]\hfill \end{array}$	$[0.001473,-1.00011]$
13	$[0.001473,-1.00011]$	$\begin{array}{c}{x}_1\in [-0.098530.101473]\hfill \\ x_2\in [-1.1001-0.90011]\hfill \end{array}$	$x_1=0,x_2=-1.00001$	N/A

GP: Goldstein and Price.

As can be seen from Tables 1 and 2, the second-level optimization is conducted from the 10th iteration. Carry out variable space updating for the optimal value obtained in iteration 9 according to formula (4) and formula (5), then take ASA algorithm to search for the optimal value in this space based on the three approximate models, obtain the updated variable space according to step 7, and then obtain the updated optimal value after further searching optimization. Iteration 10 does not meet the termination criterion yet, so continue searching in accordance with the above steps. Until termination criterion is met in the 12th iteration, the search process terminates. Take the optimal value, for which ASA algorithm searches based on QRSM model, as the global optimum. The optimum and the updated variable space in the second-level optimization are seen in Figure 3.

Figure 3.

The illustration of the second optimization stage.

Because the updated space is too small in some iterations, Figure 3 just shows the expanded and the updated variable space after searching optimization for iteration 9 and iteration 11 for convenience.

It is easy to find that the optimal point $x_1=-0.144050,x_2=-0.807620$ in the first-level optimization is close to one of the four local minimum points $x_1=-0.4,x_2=-0.6$ for GP function. However, it avoids converging to the local minima and finds the global optimum of GP function after the second-level optimization. Therefore, the proposed hybrid metamodel-based two-level global optimization method can escape the local minima and find the global optimum successfully.

Test and verify of HMTGO

Function tests

In order to verify the validity of the proposed method, this article randomly selects three two-dimensional (2D) functions, one six-dimensional (6D), and one 16-dimensional (16D) benchmark functions to test the searching ability of this method in low-dimensional optimization. Make a comparison with existing global optimization methods mode-pursuing sampling (MPS), ³⁴ HAM, ²⁵ and SAO ⁵¹

• 1. Six-hump camel-back (SC) function

$$\begin{array}{l}f(x_1,x_2)=4{x_1}^2-2.1{x_1}^4+\frac{1}{3}{x_1}^6+x_1{x}_2\\ -4{x_2}^2+4{x_2}^4{x}_1,x_2\in \left[-2,2\right]\end{array}$$

• 2. GP function

$$\begin{array}{l}f(x_1,x_2)=(1+{(x_1+x_2+1)}^2\\ (19-14x_1+3{x_1}^2-14x_2+6x_1{x}_2+3{x_2}^2))\\ ({(30+2x_1-3x_2)}^2\\ (18-32x_1+12{x_1}^2+48x_2-36x_1{x}_2+27{x_2}^2))\\ x_1,x_2\in \left[-2,2\right]\end{array}$$

• 3. Himmelblau (HM) function

$$\begin{array}{l}f(x_1,x_2)={({x_1}^2+x_2-11)}^2\\ +{(x_1+{x_2}^2-7)}^2{x}_1,x_2\in \left[-6,6\right]\end{array}$$

• 4. Hartman (HN) function

$$\begin{array}{c}f\left(x\right)=-\sum _{i=1}^4{c}_i\exp [-\sum _{j=1}^6{a}_{\mathrm{ij}}{(x_j-p_{\mathrm{ij}})}^2]x_i\in [0,1],\hfill \\ i=1,\dots ,6\hfill \\ \left[a_{\mathrm{ij}}\right]=\left[\begin{array}{cccccc}10\hfill & 3\hfill & 17\hfill & 3.5\hfill & 1.7\hfill & 8\hfill \\ \multicolumn{1}{c}{0.05}& 10\hfill & 17\hfill & 0.1\hfill & 8\hfill & 14\hfill \\ \multicolumn{1}{c}{3}& 3.5\hfill & 1.7\hfill & 10\hfill & 17\hfill & 8\hfill \\ \multicolumn{1}{c}{17}& 8\hfill & 0.05\hfill & 10\hfill & 0.1\hfill & 14\hfill \end{array}\right]\hfill \\ c_i=\left[\begin{array}{cccc}1\hfill & 1.2\hfill & 3\hfill & 3.2\hfill \end{array}\right]\hfill \\ \left[p_{\mathrm{ij}}\right]=\left[\begin{array}{cccccc}0.1312\hfill & 0.1696\hfill & 0.5569\hfill & 0.0124\hfill & 0.8283\hfill & 0.5886\hfill \\ \multicolumn{1}{c}{0.2329}& 0.4135\hfill & 0.8307\hfill & 0.3736\hfill & 0.1004\hfill & 0.9991\hfill \\ \multicolumn{1}{c}{0.2348}& 0.1451\hfill & 0.3522\hfill & 0.2883\hfill & 0.3047\hfill & 0.6650\hfill \\ \multicolumn{1}{c}{0.4047}& 0.8828\hfill & 0.8732\hfill & 0.5743\hfill & 0.1091\hfill & 0.0381\hfill \end{array}\right]\hfill \end{array}$$

• 5. 16 variables function (F16)

$$\begin{array}{c}f\left(x\right)=\sum _{i=1}^{16}\sum _{j=1}^{16}{a}_{\mathrm{ij}}(x_i^2+x_i+1)(x_j^2+x_j+1)x_i,x_j\in [-1,1],\hfill \\ i,j=1,2,\dots ,16\hfill \\ \left[a_{\mathrm{ij}}\right]=\left[\begin{array}{cccccccccccccccc}1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]\hfill \end{array}$$

Take the calculation of GP, HM, and HN as examples to study the influence of different $\rho$ for the optimization results in step 3. The searching optimization time mainly depends on the number of functions evaluations (NFE); however, Gu et al. ²⁵ claimed that when function evaluation and searching iteration are decoupled, if parallel computation is available, the number of iterations (NIT) can better represent the search efficiency. Therefore, this article adopts both the NFE and the NIT to measure the optimization efficiency of these methods. The minimum (f_min) is used to test the optimization accuracy. Table 3 lists the optimization results of different $\rho$ .

Table 3.

Comparison of the results for different $\rho$ .

Function	n	Analytical minimum	NIT			NFE			fmin
			40%	50%	60%	40%	50%	60%	40%	50%	60%
GP	2	3	9	11.6	13	50.4	63.8	68.2	19.2	3	3
HM	2	0	15	14	20	70.5	70	81	0	0	0.02
HN	6	−3.3224	11.2	11.6	15	272	273.2	325	3.1830	3.3198	3.3215

NIT: number of iterations; NFE: number of functions evaluations; GP: Goldstein and Price; HM: Himmelblau; HN: Hartman.

As Table 3 shows, compared to $\rho =50\%$ , although the NIT and NFE are decreased with smaller value of $\rho$ , the optimization accuracy has also been reduced. Moreover, the global optimum is easy to be missed, like GP function in Table 3. As for higher value of $\rho$ , more iterations and function evaluations are needed, yet the accuracy is only improved slightly or even reduced. Therefore, 50% is chosen for $\rho$ in this article.

For easy to compare, 10 runs on each of these benchmark functions are carried out in this article like MPS, HAM, and SAO. By comparison, Table 4 lists the optimization results obtained from the proposed method in this article and the existing methods.

Table 4.

Comparing of test results for low-dimensional functions.

Function	n	Analytical minimum	NIT				NFE				fmin
			MPS a	HAM a	SAO b	HMTGO	MPS	HAM	SAO	HMTGO	MPS	HAM	SAO	HMTGO
SC	2	−1.0316	11.7	6.7	–	8.7	41.2	49.6	35.7	46	−1.0308	−1.0310	−1.0315	−1.0316
GP	2	3	53.8	23.4	–	11.6	188.9	165.1	117.4	63.8	23.9077	3.0006	3.0011	3
HM	2	0	16.1	11.3	–	14	57	81.2	43.5	70	0.0388	0.0053	0.0013	0
HN	6	−3.3224	100.8	20.7	–	11.6	929.9	132	113.9	273.2	−3.2878	−3.3069	−3.3078	−3.3198
F16	16	25.8750	131.8	43.2	–	14.5	3298.2	250.6	166.8	1930.5	33.8085	26.3363	25.8970	25.987

NIT: number of iterations; NFE: number of functions evaluations; HMTGO: hybrid metamodel–based two-level global optimization method; SC: six-hump camel-back; GP: Goldstein and Price; HM: Himmelblau; HN: Hartman; F16: 16 variables function.

a

The data previously presented in the literature. ²⁵

b

The data previously presented in the literature. ⁵¹

As shown by the results in Table 4, for HM, GP, and SC, the proposed method is better than MPS on accuracy and NIT; compared with HAM method, HMTGO has advantages in accuracy and the NFE, especially for GP function, the proposed method has obvious advantages. In comparison with SAO method, although HMTGO needs more numbers of functions evaluations for SC function and HM function, the optimization accuracy is improved. In addition, the NFE and the optimization accuracy still have obvious advantages for GP function. For HN function and F16 function, the proposed method is better than MPS and HAM on optimization accuracy. Compared with SAO, the optimization accuracy is roughly equal. Meanwhile, HMTGO needs more numbers of functions evaluations than HAM and SAO, but less than MPS.

Next, one 20-dimensional (20D) and one 30-dimensional (30D) Rastrigin function are used to test the searching ability of HMTGO in high-dimensional problems. The results are presented in Table 5. Rastrigin function

$$\begin{array}{l}f(x)=10d+{\displaystyle \sum _{i=1}^d[{x_i}^2-10\cos (2\pi x_i)]}\\ x_i\in \left[-5.12,5.12\right],d=20,30\end{array}$$

Table 5.

Comparing of test results for high-dimensional functions.

Function	n	Analytical minimum	NIT	NFE	fmin
Rastrigin function	20	0	15	3151.5	0.0694
	30	0	15.7	7196.3	0.0221

NIT: number of iterations; NFE: number of functions evaluations.

As seen in Table 5, the proposed method finds the optimal value with a few iterations. However, it is noticed that the NFE non-linearly grows with the increase in the spatial dimension in Tables 4 and 5. It restricts the optimization efficiency of the proposed method in some extent. As is stated above, if parallel computation is available, function evaluation and searching iteration are decoupled, and the NIT can better represent the search efficiency. ²⁵ Meanwhile, the influence of the NIT and the NFE on the search efficiency is related to the amount of simultaneous computation. Assuming the availability of computers, the total time needed by the optimization will be represented by the NIT, rather than the number of total functions evaluations. ³⁴ Although the proposed method needs more evaluation numbers for high-dimensional functions, the NIT for searching optimization is greatly decreased, so the optimization efficiency is also improved for high-dimensional functions.

As shown in the test results of several benchmark functions, the proposed method is suitable for multi-dimensional problems. It successfully finds the global optimum in 10 trials and avoids being trapped into the local minima. Compared with the existing global optimization methods, the proposed method improves the optimization accuracy with the reduction of the NIT, which also maintains a balance between the optimization accuracy and efficiency. The test results demonstrate the superior performance of the proposed algorithm over existing search methods.

Different from the method that improves the model accuracy through sequential sampling, HMTGO method needs to reconstruct approximate model in each updated variable space. It is not difficult to find from formula (1) that the sample points needed to construct approximate models grow non-linearly with the increase in the number of design variables. Especially for the searching optimization of high-dimensional design problems, a large number of new sample points are needed to reconstruct approximate models in each iteration, which leads to a small NIT will need a large number of function evaluations.

Optimum design of tension/compression spring

The optimum design of the tension/compression spring is adopted to explain the practical application of HMTGO method in engineering field.

The popular test problem is proposed by Arora. ⁵² As shown in Figure 4, many researchers have used it as a benchmark problem in the structural optimization.^11,53,54 The minimization of the weight of a tension/compression spring is considered subject to constraints on minimum deflection, shear stress, surge frequency, and limits on the outside diameter and design variables. The design variables $x_1$ , $x_2$ , and $x_3$ are the wire diameter, the mean coil diameter, and the number of active coils.

Figure 4.

The illustration of tension/compression spring.

The mathematical model is as follows

$$\begin{array}{c}\mathrm{Min}f\left(x\right)=(2+x_3)x_1^2{x}_2\hfill \\ \mathrm{s}.\mathrm{t}.g_1\left(x\right)=1-\frac{x_2^3{x}_3}{71785x_1^4}\le 0\hfill \\ g_2\left(x\right)=\frac{4x_2^2-x_1{x}_2}{12566(x_2{x}_1^3-x_1^4)}+\frac{1}{5108x_1^2}-1\le 0\hfill \\ g_3\left(x\right)=1-\frac{140.45x_1}{x_2^2{x}_3}\le 0\hfill \\ g_4\left(x\right)=\frac{x_1+x_2}{1.5}-1\le 0\hfill \end{array}$$

where $x_1\in [0.05,2],x_2\in [0.25,1.3],\mathrm{and}{x}_3\in [2,15]$ .

It can be clearly observed from the mathematical model that the target and constraints are highly non-linear, so it is considered as a time-consuming problem. Construct QRSM, KRG, and RBF approximate models for $f(x)$ , $g_1(x)$ , $g_2(x)$ , $g_3(x)$ , and $g_4(x)$ , respectively. For the convenience of comparison, this article also carries out 11 trials for this problem. Table 6 lists the average NFE, the average of the minimum value, the minimum value of 11 trials, and the corresponding value of design variables and constraints.

Table 6.

Comparison of the results for the practical example.

Design variable	CAC Coello a	T Ray and P Saini a	S Kitayama et al. a	HMTGO
$x_1$	0.051480	0.050417	0.050000	0.060133
$x_2$	0.351661	0.321532	0.314777	0.564400
$x_3$	11.632201	13.979915	14.650042	5.443000
$g_1\left(x\right)$	−0.002080	−0.001926	−0.018820	−0.042590
$g_2\left(x\right)$	−0.000110	−0.012944	−0.006566	−0.045710
$g_3\left(x\right)$	−4.026318	−3.899430	−3.837790	−7.630490
$g_4\left(x\right)$	−0.731239	−0.752034	−0.756815	−0.583640
$f{\left(x\right)}_{11-\min }$ b	0.012705	0.013060	0.013103	0.015190
$f{\left(x\right)}_{11-\mathrm{avg}}$ c	0.012769	0.013436	0.013273	0.015248
NFE	900,000	1291	66	86.5

HMTGO: hybrid metamodel–based two-level global optimization method; NFE: number of functions evaluations.

a

The data previously presented in the literature. ¹¹

b

The minimum value of 11 trials.

c

The mean of 11 trials.

The results can be seen from Table 6 that the proposed method can obtain the optimum value which satisfies the constraints with greatly decreasing the NFE. Meanwhile, it is not difficult to find that the optimization accuracy of HMTGO is worse than the existing methods. It depends on when the constraints are considered in some extent. The existing literature ³² pointed out that if constrains are explicitly given or are of first-order or second-order for which accurate surrogates are easy to obtain, the optimum found from the approximate models of target and constraints is guaranteed to satisfy the constraints which are not approximated. For highly non-linear and computation-intensive constraints, if the space reduction strategy is adopted, then whether the real constraints are satisfied will depend on the accuracy of approximate models for constraints in the reduced variable space. In the process of reducing variable space, when the approximate models of constraints can meet the accuracy requirement, when the constraints are considered with the decrease in the variable space, and whether the real constraints are still satisfied when the variable space is small enough become problems which will have an impact on the final optimization results. This article adopts space reduction strategy to reduce the variable space and then selects partial sample points with the lowest response values under the premise of satisfying all constraints. The optimization result shows that each corresponding design variable of the found optimum values meets the real constraints. However, the accuracy is slightly lower than the existing methods. Therefore, in the global optimization process based on metamodel and space reduction strategy, when the constraints should be considered or how to deal with the constraints to obtain the optimal value with higher accuracy while satisfying all the real constraints becomes an urgent problem to be solved in the future research.

Conclusion

This article proposes a two-level global optimization method based on hybrid metamodel. This method establishes the three approximate models of QRSM, KRG, and RBF for time-consuming problems. After obtaining the approximate response value of sample points, the space reduction strategy is utilized to select partial sample points with the lowest response values to reduce variable space and guide the search to the promising region. Then, ASA algorithm is integrated with the latest approximate models to update the optimization space and obtain the final optimal value. On the basis of the hybrid metamodel and two-level optimization strategy, the optimization accuracy is guaranteed; with space reduction strategy based on hybrid metamodel, the search process is simplified and the optimization efficiency is improved. Several benchmark functions are adopted to test the searching ability of this method, and the optimum design of the tension/compression spring explains the effectiveness of the proposed method in engineering application.

As shown in the experimental results, the NFE grows exponentially with the increase in function dimension. It is mainly influenced by the minimum sample size of constructing QRSM model. Therefore, how to decrease the NFE by selecting approximate models reasonably or reducing the NIT remains to be solved. For optimization problems subject to constraints, how to deal with the constraint conditions in the process of optimization still needs further study, which will be the direction of further researches.