Surrogate-Based Optimization

Surrogate-based optimization (Queipo et al. 2005, Simpson et al. 2008) represents a class of optimization methodologies that make use of surrogate modeling techniques to quickly find the local or global optima. It provides us a novel optimization framework in which the conventional optimization algorithms, e.g. gradient-based or evolutionary algorithms are used for sub-optimization(s). Surrogate modeling techniques are of particular interest for engineering design when high-fidelity, thus expensive analysis codes (e.g. Computation Fluid Dynamics (CFD) or Computational Structural Dynamics (CSD)) are used. They can be used to greatly improve the design efficiency and be very helpful in finding global optima, filtering numerical noise, realizing parallel design optimization and integrating simulation codes of different disciplines into a process chain. Here the term “surrogate model” has the same meaning as “response surface model”, “metamodel”, “approximation model”, “emulator” etc. This chapter aims to give an overview of existing surrogate modeling techniques and issues about how to use them for optimization.


Introduction
Surrogate-based optimization (Queipo et al. 2005, Simpson et al. 2008) represents a class of optimization methodologies that make use of surrogate modeling techniques to quickly find the local or global optima. It provides us a novel optimization framework in which the conventional optimization algorithms, e.g. gradient-based or evolutionary algorithms are used for sub-optimization(s). Surrogate modeling techniques are of particular interest for engineering design when high-fidelity, thus expensive analysis codes (e.g. Computation Fluid Dynamics (CFD) or Computational Structural Dynamics (CSD)) are used. They can be used to greatly improve the design efficiency and be very helpful in finding global optima, filtering numerical noise, realizing parallel design optimization and integrating simulation codes of different disciplines into a process chain. Here the term "surrogate model" has the same meaning as "response surface model", "metamodel", "approximation model", "emulator" etc. This chapter aims to give an overview of existing surrogate modeling techniques and issues about how to use them for optimization.

Overview of surrogate modeling techniques
For optimization problems, surrogate models can be regarded as approximation models for the cost function (s) and state function (s), which are built from sampled data obtained by randomly probing the design space (called sampling via Design of Experiment (DoE)). Once the surrogate models are built, an optimization algorithm such as Genetic Algorithms (GA) can be used to search the new design (based on the surrogate models) that is most likely to be the optimum. Since the prediction with a surrogate model is generally much more efficient than that with a numerical analysis code, the computational cost associated with the search based on the surrogate models is generally negligible.
Surrogate modeling is referred to as a technique that makes use of the sampled data (observed by running the computer code) to build surrogate models, which are sufficient to predict the output of an expensive computer code at untried points in the design space. Thus, how to choose sample points, how to build surrogate models, and how to evaluate the accuracy of surrogate models are key issues for surrogate modeling.

Design of experiments
To build a surrogate model, DoE methods are usually used to determine the locations of sample points in the design space. DoE is a procedure with the general goal of maximizing the amount of information gained form a limited number of sample points (Giunta et al., 2001). Currently, there are different DoE methods which can be classified into two categories: "classic" DoE methods and "modern" DoE methods. The classic DoE methods, such as full-factorial design, central composite design (CCD), Box-Behnken and D-Optimal Design (DOD), were developed for the arrangement of laboratory experiments, with the consideration of reducing the effect of random error. In contrast, the modern DoE methods such as Latin Hypercube Sampling (LHS), Orthogonal Array Design (OAD) and Uniform Design (UD) (Fang et al., 2000) were developed for deterministic computer experiments without the random error as arises in laboratory experiments. An overview of the classic and modern DoE methods was presented by Giunta et al. (2001). A more detailed description of existing DoE methods is beyond the scope of this chapter.
The schematics of 40 sample points selected by LHS and UD for a two-dimensional problem are sketched in Figure 1.

Surrogate models
There are a number of surrogate models available in the literatures. Here we limit our discussion to three popular techniques such as RSM (polynomial Response Surface Model), Kriging, RBFs (Radial Basis Functions).
For an m-dimensional problem, suppose we are concerned with the prediction of the output of a high-fidelity, thus expensive computer code, which is correspondent to an unknown function m : The pair ( S , S y ) denotes the sampled data sets in the vector space.
With the above descriptions and assumptions, our objective here is to build a surrogate model for predicting the output of the computer code for any untried site x (that is, to estimate () y x ) based on the sampled date sets ( S , S y ), in an attempt to achieve the desired accuracy with the least possible number of sample points.

Quadratic response surface method
Here we use "RSM" to denote a polynomial approximation model in which the sampled data is fitted by a least-square regression technique. In RSM-based optimization applications, the "quadratic" polynomial model usually provides the best compromise between the modeling accuracy and computational expense, when compared with the linear or higher order polynomial models. An advantage of RSM is that it can smooth out the various scales of numerical noise in the data while captures the global trend of the variation, which makes it very robust and thus well suited for optimization problems in engineering design.
The true quadratic RSM can be written in the following form: where () y x is the quadratic polynomial approximation and ε is the random error which is assumed to be normally distributed with mean zero and variance of 2 σ . The error i ε at each observation is supposed to be independent and identically distributed. The quadratic RSM predictor () y x can be defined as: where 0 β , i β , ii β and ij β are unknown coefficients to be determined. Since there are totally (1 ) (2 ) / 2 pm m =+ + unknown coefficients in Eq.(4), building a quadratic RSM with m variables requires at least p sample points. Let p ∈ β  be the column vector contains these p unknown coefficients. The least square estimator of β is where After the unknown coefficients in β are determined, the approximated response ŷ at any untried x can be efficiently predicted by Eq. (4).

Kriging model
Different from RSM, Kriging (Krige, 1951) is an interpolating method which features the observed data at all sample points. Kriging provides a statistic prediction of an unknown function by minimizing its Mean Squared Error (MSE). It can be equivalent to any order of polynomials and is thus well suited for a highly-nonlinear function with multi extremes. For the derivation of Kriging (Sacks et al., 1989), the output of a deterministic computer experiment is treated as a realization of a random function (or stochastic process), which is defined as the sum of a global trend function T () fx β and a Gaussian random function () Z x as following where T 01  denotes the vector of the corresponding coefficients. In general, T () fx β is taken as either a constant or low-order polynomials. Practice suggests that the constant trend function is sufficient for most of the problems. Thus, T () fx β is taken as a constant 0 β in the text hereafter. In Eq. (7), () Z ⋅ denotes a stationary random process with zero mean, variance 2 σ and nonzero covariance of Here (, ) R ′ xx is the correlation function which is only dependent on the Euclidean distance between any two sites x and ′ x in the design space. In this study, a Gaussian exponential correlation function is adopted, and it is of the form  From the derivation by Sacks et al. (1989) the Kriging predictor () y x for any untried x can be written as T1 00 where the generalized least square estimation of 0 and n ∈ 1  is a vector filled with ones, and , Rr are the correlation matrix and the correlation vector, respectively. and Rr are defined as A unique feature of Kriging model is that it provides an uncertainty estimation (or MSE) for the prediction, which is very useful for sample-points refinement. It is of the form Assuming that the sampled data are distributed according to a Gaussian process, the responses at sampling sites are considered to be correlated random functions with the corresponding likelihood function given by The optimal estimation of 0 β and the process variance are obtained analytically, yet depends on the remaining hyper-parameters which can be solved by a numerical optimization algorithm such as GA. The hyperparamters tuning strategies are discussed by Toal et al. (2008). Note that the above Kriging formulation can be extended by including gradient information obtained by Adjoint method (Han et al. 2009, Laurenceau et al. 2008 or lower-fidelity data by lower-fidelity analysis code , Forrester et al. 2007).

Radial basis functions
In additional to Kriging, RBFs model (Hardy, 1971) is known as an alternative interpolation method for surrogate modeling. For the RBFs approach by Powell (1987), the approximation of the unknown function () y x at an untried x is formally defined as the linear combination of the radial basis functions and a global trend function as where i ω are the i-th unknown weight coefficient, x are the basis functions that depend on the Euclidean distance between the observed point () i x and the untried point x (similar to the correlation function of kriging model); () P x is the global trend function which is taken as a constant 0 β here. To ensure the function values at observed points are reproduced by the RBFs predictor, the flowing constraints should be satisfied: Then the additional constraints for () P x should be imposed as Solving the linear equations formed by Eq. (18) Where T1 1 T1 When the above RBFs predictor is compared with the Kriging predictor (see Eq. (10)), one can observe that they are essentially similar, only with the basis-function matrix Ψ (also called Gram matrix) and the basis function vector () φ x being different from the correlation matrix R and the correlation vector () rx of the Kriging predictor, respectively. In addition, RBFs differs from Kriging at the following two aspects: 1) RBFs doesn't provide the uncertainty estimation of the prediction; 2) The model parameters can't be tuned by MLE like Kriging. Generally, Kriging can be regarded as a particular form of RBFs.
To build a RBFs model, one needs to prescribe the type of basis functions that only depends on the Euclidean distance r ′ =− xx between any two sites x and ′ x . Compared to the correlation function used for a Kriging model, more choices are available for a RBFs model, which are partially listed in Table 1.  Table 1 can be classified into two categories: decaying functions (such as GAUSS and HIMQ) and growing functions (POW, TPS and HMQ). The decaying functions can yield positive definite matrix Ψ , which allows for the use of Cholesky decomposition for its inversion; the growing functions generally result in a nonpositive definite matrix Ψ and thus LU decomposition is usually used alternatively. The schematics of the basis functions for one-dimensional problem is sketched in Figure 3.

The big (ger) picture
In addition to what we mentioned above, there are also a few surrogate models available in the literatures, such as Artificial Neutral Network (ANN) (Elanayar et al. 1994, Park et al. 1991, Multiple Adaptive Regression Splines (MARS) and Support Vector Regression (SVR) (Smola & Schoelkopf 2004). Although these methods are coming from different research communities, the idea is similar when using them for function prediction in surrogate modeling. They are not described in detail here due to the limited space. The readers are referred to read the paper by Wang & Shan (2007) and the books written by Keane et al. (2005) and by Forrester et al. (2008) for more description about surrogate modeling techniques.

Evaluation of approximation models
An important issue for the surrogate modeling is how to estimate the error of the approximation models. Only when the surrogate model with sufficient accuracy is built can the reliable optimum be obtained. Here we use two variables ( e and e σ ) to evaluate the error of the approximation models at test points, which are also chosen by DoE method. The average relative error is

Framework of building surrogate models
A Generic framework of building a surrogate model is sketched in Figure 4. Note that the initial surrogate model can be evaluated by Eqs. (22) and (23) and a branch for resampling is denoted by black dashed line in Figure 4.

Use of surrogate models for optimization
Suppose we are concerned with solving a general optimization problem as Objective minimize y() where c n is the number of state functions which is in line with the number of inequality constraints (assuming that all the equality constraints have been transformed into inequality constraints.); l x and u x are the lower and upper bound of design variables, respectively; the object function y() x and state functions g () i x are evaluated by an expensive analysis code. Traditionally, the optimization problem is solved by either a gradient-based algorithm or a gradient-free algorithm such as GA. It may become prohibitive due to the large computational cost associated running the expensive analysis code. Alternatively, here we are concerned with using surrogate modeling techniques to solve the optimization problem, in an attempt to dramatically improve the efficiency.

A simple framework
The basic idea of using surrogate models in optimization can be quite simple. First, the surrogate models for the object function(s) and state function(s) with sufficient accuracy are built (see Figure 2); second, the optimum is found by an optimizer, with the object function(s) or state function(s) evaluated by surrogate models, rather than by the expensive analysis code. Since prediction with the surrogate models is much more efficient than that by the expensive analysis code, the optimization efficiency can be greatly improved. The comparison of the conventional optimization and surrogate-based optimization is sketched in Figure 5. In addition to improve optimization efficiency, surrogate models also serve as an interface between the analysis code and the optimizer, which makes the establishment of an optimization procedure much easier. One of the examples for such a surrogate-based optimization framework can be found in paper by .

A bi-level framework
Although the framework of the surrogate-based optimization sketched in Figure 5. (b) is very intuitive and simple, questions may arise: are the surrogate models accurate enough? has the true optimum been reached? In fact, the optimum gained by the surrogate models is only an approximation to the true optimum (see Figure 5. (a)). One has to refine the surrogate models by adding new sample points, which is to be observed by running the analysis code. The procedure of augmenting new sample point(s) to the current sampled data sets is the so-called "sample-point refinement". The rules of determining the new sample sites towards the true optimum are called "infill criteria", which will be discussed in section 3.2. The flowchart of a surrogate-based optimization with additional process of sample-point refinement is sketched in Figure 6. It can be regarded as a bi-level optimization framework, with the process of building and refining the surrogate models (which needs to run the expensive analysis code) acting as the main optimization and the process of using surrogate models to determine the new sample sites acting as the sub-optimization(s).

Infill criteria
The infill criterion is related to the determination of the new sample sites by solving suboptimization problem(s). Three infill criteria are discussed here: Searching Surrogate Model(s) (SSM), Expected Improvement (EI) and Lower-Bounding Confidence (LCB).

Searching surrogate model(s)
Provided that the initial surrogate models have been built, an optimizer such as GA can be used to find the optimum, which in turn can be employed to refine the surrogate models.
The mathematical model of the sub-optimization for determining the new sample site is of the form Ôbjective minimize y() where ŷ() x and ĝ () i x are surrogate models of y() x and g () i x , respectively. With the optimal design variables opt x gained by the surrogate models in hand, one needs to run the expensive analysis code to compute the corresponding true function value and compare it with what predicted by the surrogate models. If the error between them is blow a threshold, the optimization process can be terminated; if not, the new sample point is augmented to the sampled data sets and the surrogate models are rebuilt; the process is repeated until the optimum solution is approached.
This criterion applies for all the surrogate models and is very efficient for local exploitation of the promising region in the design space.

Expected improvement
Surrogate model such as Kriging provides not only an function prediction but also an estimation of the mean squared error (MSE). In fact, the prediction by a Kriging model, ˆ( ) y x , at any point can be regarded as a Gaussian random variable with the mean given by the Kriging predictor, and the variance given by the mean squared error, () 2 s x (see section 2.2.2). Viewed in this way, a probability can be computed that the function value at any untried x would fall below the minimum among the sample points observed so far. Then Expected Improvement (EI) function (Jones et al 1998, Jeong et al. 2005 can be calculated to account for the improvement of the object function we expect to achieve at any untried x . The definition of EI is of the form where () Φ  and () φ  are the cumulative distribution function and probability density function of a standard normal distribution, respectively.

Real-World Applications of Genetic Algorithms 354
denotes the minimum of the observed data so far. The greater the EI, the more improvement we expect to achieve. The point with maximum EI is located by a global optimizer such as GA then observed by running the analysis code. For this infill criterion, the constraints can be accounted by introducing the probability that the constraints are satisfied. The corresponding sub-optimization problem can be modeled as where () x i s denotes the estimated standard error corresponding to the surrogate model The optimum site opt x obtained by solving Eq. (27) is observed by running analysis code and the new sample point is added to the sampled date sets; the surrogate models are rebuilt and the whole process is repeated until the global optimum is approached.

Lower-bounding confidence (LCB)
The LCB function is defined as the weighted sum of predicted function value ˆ( ) y x and the standard error of the prediction ˆ( ) s x . For an optimization problem of finding the minimum of the unknown function () y x , a simple expression for LCB function is of the form where A is a constant which balances the influence of the predicted function and the corresponding uncertainty. Best practice suggests 1 A = works well for a number of realistic problems. The corresponding sub-optimization problem can be modeled as The above optimization problem can be solved via a global optimizer such as GA. Since the point with smallest value of LCB indicates the possible minimum of the unknown function, the optimum site opt x is then observed and added to sampled data sets to refine the surrogate models. This procedure is performed iteratively until the global optimum is reached.

Airfoil design
Using an in-house Navier-Stokes flow solver, the objective of the problem is to minimize the drag of an RAE2822 airfoil at the flow condition of where Area 0 ,C l0 , C m0 are the area, lift coefficient, and moment coefficient of the baseline airfoil, respectively. The first constraint is in consideration of the structural design of the wing to guarantee the volume of the wing; the second one is to enforce a constant lift of the wing in order to balance the weight of the aircraft at cruise condition; the third one is to control the pitching moment of the airfoil to avoid large drag penalty of the horizontal tail paid for balancing the aircraft.
The initial number of samples for Kriging is set to 20, selected by the Latin Hypercube Sampling (LHS). The airfoil is parameterized by 10 Hicks-Henne bump functions (Hicks & Henne, 1978); and the maximum amplitude of each bump is max / 0.544% Ac = . Both of the SSM and EI infill strategies are adopted in the surrogate refinement. Table 2 presents the optimization results of the two optimization method. The optimized and initial airfoils and the corresponding pressure coefficient distributions are compared in Figure 7. Note that the aerodynamic coefficients of the initial airfoil RAE2822 are set to 100. Obviously, the Krigingbased optimization method gives better result, and with higher efficiency, and is more likely to find the global optimum.

Wing design
Here we are concerned with the preliminary design for a high-subsonic transport-aircraft wing of a wing/body combination, considering aerodynamic, structure and static aeroelastic effect. The calculation of the external flow is carried out by numerical solutions of the full potential equation in conservative form (Kovalev & Karas, 1991). The FEM-based commercial software ANSYS is used for analyzing the structural performance of the wing with double-beams sandwich structure. Weak coupling method is adopted for static aeroelastic analysis.
The optimization objectives are to maximize the aircraft lift-to-drag ratio and minimize the weight of wing for a fixed maximum take-off weight of 54 tons and cruise Mach number of 0.76 at 10,000 meters high. The wing is composed of inner and outer wing. The reference area of wing is 105 square meter. The mathematical model for optimization is of the form Eight supercritical airfoils are configured along the span. The optimization is subject to 4 constraints. The first constraint is to enforce a constant lift of the wing in order to balance the weight of the aircraft at cruise condition; the second one is to guarantee a near constant wing loading; the third and fourth constraints are to make sure that the strength and rigidity requirements are satisfied. The definition for the limits of design variables is listed in Table  3. The first four design variables define the aerodynamic configuration of the wing. The four remain are for structure design. The detail can be found in paper by Zhang et al. (2008).
The uniform design table, U100(10 8 ), is used to creates one hundred candidate wings for building surrogate models. The other forty-five candidate wings are created by the uniform design table, U45(10 8 ), for evaluating the approximation model. For each wing, the static aeroelastic analysis is performed to obtain the responses of lift (L), lift-to-drag ratio (L/D), wing area (S wing ), maximum stress (σ max ), maximum deformation (δ max ) and wing weight (W wing ). Then the average relative errors and the root mean squared errors are calculated to evaluate the approximation models, as listed in Thickness of back-spar web mm 2 6 Thickness of lower skin mm 3 7 Thickness of Upper skin mm 3 7 Then the multi-objective optimization for the supercritical wing is performed based on RSM due to its higher computational efficiency. Weighted sum method is used to transform the multi-objective optimization into a single-objective optimization. Sequential quadratic programming method is employed to solve the optimization. One of the candidate wings with better performance, are selected as the initial point for optimization. The optimal design is observed by running the analysis codes and the results are listed in Table 5. Where 0 X and 0 Y is the initial wing scheme and its response, respectively; * X and * Y is the optimal wing scheme and its actual response, respectively; Ŷ is the response at * X calculated by the approximation models. For the optimal wing scheme, the largest relative error of approximation models is no more than 3 percent. It again proves the high accuracy of the approximation models. Figure 8 shows the contour of the equivalent stress of the optimal wing. It shows that the stress is larger in the location of intersection of inner wing and outer wing due to the inflexion. Figure 9 shows pressure distribution of the optimal wing. It shows that the wing basically meets the design requirements of a supercritical wing. A little bit non-smoothness of the pressure distribution may be caused by non-uniform deformation of the skin. Figure  10 shows the convergence history of aeroelastic deformation, which shows that fast convergence of the aeroelastic deformation of the optimal wing.
The optimization , together with the aeroelastic analysis of all candidate wings, only takes about two days on a personal computer of Pentium(R) 4 CPU 2.8GHz. If more computers are used to concurrently calculate the performance of different candidate wings, the cost can be further greatly reduced.

Conclusion
An overview of the existing surrogate models and the techniques about how to use them for optimization is presented in this chapter. Among the surrogate models, the regression model such as the quadratic response surface model (RSM) is well suited for a local optimization problem with relatively simpler design space; interpolation models such as Kriging or RBFs can be used for highly-nonlinear, multi-modal functions, and thus well suited for a global problem with relatively more complicated design space. From an application point of view, the simple framework of surrogate-based optimization is a good choice for an engineer design, due to the fact that surrogate model can act as an interface between the expensive analysis code and the optimizer and one doesn't need to change the analysis code itself. The drawback of this framework is that the accuracy of optimum only depends on the approximation accuracy of surrogate model and we generally get an approximation to the true optimum. In contrast, the bi-level framework with different infill criteria provides an efficient way to quickly find true optimum without the need of building globally accurate surrogate models. Multiple infill criteria seem to be a better way to overcome the drawback of the single infill criterion.
Examples for airfoil and wing designs show that surrogate-based optimization is very promising for aerodynamic problem with number of design variables being less than about 10. For higher-dimensional problem, the computational cost increases very quickly, which can be prohibitive. Thus, use of surrogate model for high(er)-dimensional optimization problems would become an important issue of future work.