Keywords: expensive optimization problem, evolution strategy, CMA-ES, Kriging model, pre-selection, evolution control, approximate ranking procedure, parameter identification. V problè me d'optimisation coûteuse, straté gie d'é volution, CMA-ES, modè le de Krigeage, pré sé lection, contrôle d

Mayssam Al Soufi , Hao BAI, Yujun CAO, Zelong MA, ……, thank you very much for your valuable supports and aids in

XIII

List of Figures

 Table 4.
General Introduction

Optimization is extensively required and applied in almost all disciplines, whether economics, sciences, or engineering. Commonly, in engineering design and management, we almost always try to optimize something, for instance, to minimize the cost and energy consumption, or to maximize the profit, performance and efficiency. In real life, resources, time and money are always limited; consequently, optimization is far more important.

Motivated by industrial and research demands, many powerful optimization techniques (algorithms) have been developed. The algorithm chosen for an optimization task largely depend on the type of the problem, the nature of the algorithm, the desired quality of solution, the available computing resource, time limit, availability of the algorithm implementation, and the expertise of the decision-makers [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]. This thesis devotes to a class of challenging optimization problems known as expensive optimization problems [START_REF] Tenne | Computational Intelligence in Expensive Optimization Problems[END_REF]. In this introductory chapter, the scope and motivation of our study are firstly explained. Then, we give a brief outline of the whole thesis.

Aim and Scope

With the development of numerical methods and computing technology, contemporary engineering design is heavily based on computer simulations, such as finite element method (FEM) and computational fluid dynamics (CFD) simulations. Correspondingly, computeraided design optimization is now involved in a wide range of applications, for instance, optimization of automotive components and shape optimization of wind turbine blades.

However, computer-aided design optimization is accompanied by several difficulties, which bring high computational cost in optimization. Thus, computer-aided design optimization generally belongs to the so-called expensive optimization problems. The computational challenge arises in this class of expensive optimization problems generally are due to:

Aim and Scope 2  Objective function is evaluated based on computer simulation, which can be expensive (require anywhere from minutes to hours even days of computation time for each run of simulation). Therefore, in this case, the computational cost would be very expensive, because usually a large number of objective function evaluations is needed in the optimization process.

 There is no explicit analytical expression for the objective function or its derivatives. Thus, derivative-free algorithms are required to solve this black-box type of optimization problem. Generally, without the aid of derivative information, derivative-free algorithms require more evaluations of objective function than derivative-based algorithms.

 For many computer-aided design optimization problems in industry and engineering, the landscape of objective function may be non-smooth, multimodal, discontinuous and ill-conditioned. These difficulties also bring about high computational cost in finding the optimum.

In this thesis, our objective is to develop powerful optimization techniques that can efficiently deal with the expensive optimization problems subject to the box constraints. This kind of problems can be described as

Thesis Outline

This thesis is organized in six chapters. After this general introduction, Chapter 1 gives the state of the art of optimization. Some basic concepts of optimization and a brief overview of optimization algorithms are firstly presented in Chapter 1. Then, a review of surrogateassisted evolutionary computation based fitness approximation is offered.

Chapter 2 systematically presents the evolution strategies (ESs), from the main principles and evolutionary operators to strategy parameter control and algorithms. The to analyze and evaluate the performance of existed and our newly developed algorithms.

Chapter 5 provides an engineering application of the KA-CMA-ES algorithm using approximate ranking procedure in material parameter identification of an elastic-plastic damage model.

Chapter 6 gives a brief summary of the whole thesis in French.

Lastly, the conclusion and perspective part concludes the thesis and summarizes our contributions to surrogate-assisted evolution strategies for expensive optimization problems.

State of the Art of Optimization

This chapter presents a brief overview of optimization techniques and the state of the art of surrogate-assisted evolutionary optimization. Firstly, basic concepts of optimization are concisely introduced in Section 1.1. An overview of optimization algorithms is given subsequently in Section 1.2, including derivative-based and derivative-free algorithms.

Section 1.3 reviews the development of surrogate-assisted evolutionary algorithms (EAs) for optimization.

Basic Concepts of Optimization

Optimization means finding the best solution, which is specified by certain goals, such as minimizing the cost of a process or maximizing the efficiency of a system from all feasible solutions [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]. Feasible solutions are those that satisfy all the conditions or constraints in the optimization problem. Nowadays, optimization is required and applied in almost all disciplines, whether economics, science, or engineering. In this section, some basic notions of optimization are presented.

In an optimization problem, a function that mathematically represents the problem is to be minimized or maximized. The function that is being optimized is referred to as the objective function. Generally, the objective function is a quantity such as cost, profit, efficiency, size, weights, output, and so on. The variables in the objective function are denoted as the design variables. For example, for a structural optimization problem, design variables could be the dimensions of a structure or its material parameters. In some problems, design variables take on values from a discrete set, often a subset of integers, whereas in other problems, design variables can take on any real values. Problems with discrete design variables are called discrete optimization problems; similarly, problems with continuous variables are continuous optimization problems. In some optimization problems, there are some constraints on design variables, which restrict the domain from which design variables

      0, 1, , and 0, 1, , , d ij S
g i m h j p       x x x (1.2)
which is called feasible set (also called feasible region) or search domain, and any element from this set is called a feasible point. A feasible point  x is called global minimum (global optimum) solution, if     for all .

f f f S     x x x (1.3)
Conversely, it is called a local minimum (local optimum) if the above inequality holds for x only within its neighborhood [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF]. The global and local optimum are simply illustrated in Figure 1.1.

The optimization problem described in formulation (1.1) is called a constrained optimization problem since there are constraints on design variables. Optimization problem without constraints are correspondingly referred to as unconstrained optimization problems.

Unconstrained optimization problems arise directly in many practical applications. They are also the basis for constrained optimization problems in which the constraints can be replaced by a penalty term in the objective function. Constrained optimization problems come from applications in which there are explicit constraints on the variables. The constraints on the variables can vary widely from simple bounds to systems of equalities and inequalities that model complex relationships among the variables. Constrained optimization problems can be furthered classified according to the nature of the constraints (e.g., linear, nonlinear, convex) and the smoothness of the functions (e.g., differentiable or non-differentiable).

Different kinds of optimization problems usually call for distinct approaches to solve them. In optimization, after an problem is modeling as a mathematical form similar as (1.1), the next essential step is to use the proper algorithm to solve the optimization problem. A brief review of optimization algorithms is given in the following section.

Overview of Optimization Algorithms 1.2.1 Classification of Optimization Algorithms

Optimization algorithms, which try to find the minimum values of mathematical functions, are essential in optimization. There are many optimization algorithms in the literature and no single algorithm is suitable for all problems. Thus, the right choice of an optimization algorithm is crucially important in solving a given optimization problem.

Optimization algorithms can be classified in many ways, depending on the focus or the characteristics we are trying to compare [START_REF] Koziel | Computational Optimization[END_REF]. For instance, optimization algorithms can be Several typical and popular algorithms including both derivative-based and derivative-free algorithms for multi-variable/multi-dimensional unconstrained optimization problems are introduced subsequently.

Derivative-based Algorithms

Since derivative-based algorithms require the first-order and second-order derivatives of the objective function, firstly, some concepts about the derivative are given. The firstorder derivative or so-called gradient of a function  

f x at a point   T 1 ,, d xx  x is         T 12 , , , . d f f f f x x x           x x x x (1.4)
The gradient represents the slope of the tangent of the graph of the function. More precisely, the gradient points in the direction of the greatest rate of increase of the function [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]. It is obvious that   f x decreases fastest in the direction of its negative gradient. Thus, gradient information provides the search direction to locate the minimum of the function.

The Hessian matrix or Hessian   Hx , which represents the second-order partial

derivatives of   f x , is written as     2 2 2 2 1 1 2 1 2 2 2 2 2 2 1 2 2 2 2 2 2 12 . d d d d d f f f x x x x x f f f f x x x x x f f f x x x x x                                       H x x (1.5)
For unconstrained design space, the conditions for optimal solution can be formulated via the first and the second partial derivatives. The necessary conditions for  x to be a minimum of   ).

In practice, for most optimization problems in which the objective functions   f x are nonlinear, the derivatives have to be evaluated numerically. Finite difference methods, including forward difference, backward difference, and central difference methods, can be used to calculate the derivative of a function at a point. Because the central difference method for computing the derivative is more accurate than forward/backward difference methods, it is preferable in applications.

Steepest Descent Method

The steepest descent or gradient descent algorithms choose search direction as the negative gradient, i.e.,

      kk f   dx
based on the idea that f decreases fastest in the direction of this negative gradient [6]. In successive iterations, the design variables are updated as

                1 , k k k k k k k f        x x d x x (1.6)
where   k  is called step length, which is a positive scalar parameter that can be determined using the line search algorithm such as the golden section method [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]. The iterative process terminates until the convergence criteria are fulfilled. The algorithm for the steepest descent method is described in Algorithm 1.1. The steepest descent method ensures a reduction in the objective function value at every iteration. Generally, when the search point is far away from the minimum, the gradient will be higher and the function reduction will be large for an iteration. As approaching the minimum, the gradient value usually decreases, i.e., the method becomes sluggish (slow convergence) near the minimum. To overcome this, the conjugate gradient method was proposed, which selects successive descent directions in a conjugate direction to previous descent directions [7].

Algorithm 1.1 Steepest Descent 1: input: objective function   f x , initial point   0 x , tolerance  , set 0 k  . 2: repeat 3: compute the search direction       kk f   dx 4: determine the step length   k  by minimizing         k k k f   xd using line search 5: update candidate         1 k k k k    x x d 6: 1 kk  7: until     k f   x

Newton's Method

Newton's method uses the search direction, which is based on the first-order and second-order derivative information, given by

          1 k k k f       d H x x .
Correspondingly, Newton's method updates candidate solutions at each iteration via

                1 1 , k k k k k k f          x x d x H x x (1.7) where   f  x and     2 f  H x
x denotes the gradient vector and Hessian matrix of   f x , respectively. The algorithm for Newton's method is described in Algorithm 1.2.

From Equation (1.7), for updating candidate solution at each iteration, the inverse of Hessian matrix need to be computed, which can be expensive. To ease the computational cost, approximations to the Hessian and its inverse are used. This has brought about the socalled quasi-Newton methods, in which the DFP and BFGS methods are popular [6,[START_REF] Papadrakakis | Survey of quasi-Newton methods with reduced storage[END_REF].

Algorithm 1.2 Newton's method 1: input: objective function   f x , initial point   0 x , tolerance  , set 0 k  . 2: repeat 3: compute the search direction           1 k k k f       d H x x 4: update candidate       1 k k k   x x d 5: 1 kk  6: until     k f   x

Levenberg-Marquardt Method

The advantage of the steepest descent method is that it reaches closer to the minimum of the function in a few iterations even when the starting point is far away from the optimum [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]. However, the method shows sluggishness near the optimum point. On the contrary, Newton's method shows a faster convergence if the starting point is close to the minimum point. Newton's method may not converge if the starting point is far away from the optimum point. These have driven the development of the so-called hybrid methods that take advantages of different methods.

Algorithm 1.3 Levenberg-Marquardt (LM) method 1: input: objective function   f x , initial point   0 x , tolerance  , set 0 k  . 2: repeat 3: compute         , kk ff  xx and     k Hx 4: compute the search direction           1 k k k f         d H x I x 5: update candidate       1 k k k   x x d 6: If        
    k f   x
The Levenberg-Marquardt (LM) method is a hybrid method that combines the strength of both the steepest descent and Newton's methods. The search direction in LM method is:

1.2 Overview of Optimization Algorithms 12           1 , k k k f         d H x I x (1.8)
where Ι is an identity matrix and  is a scalar that is set to a high value at the start of the algorithm. The value of  is altered during every iteration depending on whether the function value is decreasing or not. If the function value decreases in the iteration,  decreases by a factor. On the other hand, if the function value is increase in the iteration,  increases by a

factor. The algorithm for the LM method is described in Algorithm 1.3.

Derivative-free Algorithms

Algorithms using derivative information are generally efficient, but these algorithms pose limitations on objective functions. In many real-world optimization problems, the objective function is evaluated by performing computer experiments and there is no analytical expression for the objective function or its derivatives. In these case, derivativefree algorithms are required. Additionally, when discontinuity exists in objective function, derivative-free algorithms may be more efficient and proper. Subsequently, several derivative-free algorithms are presented.

Nelder-Mead Simplex Algorithm

The Nelder-Mead algorithm [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] solves the optimization problem by containing the solution within a simplex. Simplex refers to a geometric figure formed by 1 d  points in a ddimension space, which is the generalization of a polygon to d-dimension. For simplicity, simplex in the d-dimension space is referred to as d-simplex, which has 1 d  vertices. The Nelder-Mead algorithm starts with a set of points in d forming a simplex and at each iteration, the objective function is evaluated at the vertice of the simplex. Using this information, the simplex is moved in the search space. The process of moving the simplex is continued until the optimum value of the function is reached. Three basic operations are required to move the simplex in the search space: reflection, contraction, and expansion [START_REF] Arora | Optimization: Algorithms and Applications[END_REF].

: objective function   f x , initial simplex       1 2 1 , , , d  x x x
, coefficients ,,    and  , tolerance  .

2: repeat

3: re-order the vertices according to

            1 2 1 d f f f     x x x 4: compute the centroid point   1 1 d i d i   xx 5:
generate a trial point through reflection

      1 rd      x x x x 6: if         1 r ff  xx then 7: compute expansion point         e r r     x x x x 8: if         er ff  xx then 9: accepted the expansion point     1 de   xx 10: else 11: accept the reflection point     1 dr   xx 12: end if 13: else if         rd ff  xx then 14: accept reflection point     1 dr   xx 15: else if         rd ff  xx then 16: contraction         11 c d d      x x x x 17: if         1 cd ff   xx then 18: accept the contraction point     1 dc   xx 19: else 20: reduction (shrink)           11 , 2,3, , ii id      x x x x 21: end if 22: end if 23: until       2 1 1 1 1 d i d i ff        xx
The fundamental procedure of Nelder-Mead algorithm is described in the following.

The first step is to construct an initial d-simplex with 1 d  vertices and to evaluate the objective function at the vertices. Then, the 1 d  vertices are ordered according to their corresponding objective function values as

            1 2 1 , d f f f     x x x (1.9)
where

      1 2 1 , , , d  x x
x are ordered vertices of the simplex, whose objective function values are in ascending order. Apparently,  

1 d  x
is the worse point (solution) and   1

x is the best solution, which is the convention in simplex method for minimization problem. At each iteration, similar ranking manipulations are performed.

Then, the centroid point x of simplex is computed using all the vertices but with the exclusion of the worst vertex  

1 d  x . That is   1 1
.

d i i d    xx (1.10)
Using the centroid point as the basis point, the reflection point of the worse point  

1 d  x is computed as       1 , 0, rd       x x x x (1.11)
where  is the so-called reflection coefficient. The typically value of 1   is often used.

Whether the new trial solution is accepted or not and how to update the new vertex, depends on the objective function value at   r x . There are three possibilities [10]:

 If             1 rd f f f  x x x
, then replace the worst vertex  

1 d  x by   r x , i.e.,     1 dr   xx .  If         1 r
ff  xx which means the objective improves, then it is possible to move the vertex further along the line of reflection to seek an expand point that can improve the objective even further. The expansion point is computed as

        , e r r     x x x x (1.12)
where  is the expansion coefficient. Typically, 2   is adopted. Now, we have to

check if     e f x improves even better. If         er ff  xx , it is accepted and the vertex is updated as     1 de   xx ; otherwise, the result of reflection is used, i.e.,     1 dr   xx .
 If the reflect point is not better than the second worst point, i.e.,

        rd ff 
xx , the contraction operation, which reduce the size of the simplex while maintaining the best sides, is performed. The contraction point is computed as

        11 , 0 1, c d d        x x x x (1.13)
where  is the so-called contraction coefficient and

12   is usually used. If         1 cd ff   xx is true, we then update     1 dc   xx .
If all the above steps fail, we should reduce (shrink) the size of the simplex towards the best vertice   1 x . This is called the reduction (or shrink) step, which is expressed by

          11 , 2,3, , , ii id      x x x x (1.14)
where is called the shrink coefficient, and usually

12

  is used. The preceding operations are continued until the the standard deviation of the objective function values computed at the vertices of the simplex becomes less than the tolerance  , i.e.,

      1 2 1 1 . 1 d i i ff d         xx (1.15)
Above described Nelder-Mead simplex algorithm is described in Algorithm 1.4.

Trust-Region Method

The so-called trust-region method is one of the most widely used optimization algorithms. This method uses a model, which is usually smooth, easy to evaluate, and presumed to be accurate in a neighborhood (trust-region) about the current iteration, to approximate the objective function [START_REF] Rios | Derivative-free optimization: A review of algorithms and comparison of software implementations[END_REF]. The next trial solution is then found by using the approximate model of objective function. Since the approximate model of objective function is smooth and easy to evaluate, finding the next solution based on approximation is cheaper than that from original function.

The fundamental step in trust-region algorithm is to approximate the objective function.

The commonly used approximate model is the quadratic model [START_REF] Conn A. R | Recent progress in unconstrained nonlinear optimization without derivatives[END_REF], which is of the form:

                        T T 1 , 2 k k k k k k k f          x x g x x x x H x x (1.16)
and this approximation is valid only in a small neighborhood of   k x , which is the so-called trust-region, defined by

        , k k k d       x x x (1.



x to interpolate a set of sample points as in [START_REF] Conn A. R | Recent progress in unconstrained nonlinear optimization without derivatives[END_REF].

Let    x denotes the minimum of approximation model   k  in the trust-region   k  , that is                             T T arg min 1 . 2 k k k k k k k k f             x xx x g x x x x H x x (1.18)
How good the approximate model

    k 
x is to the actual objective function   f x in the trust-region can be measured by the ratio between the actual reduction and the predicted reduction

                      . k k k k k ff r       xx xx (1.19)
Now, the two key questions are how to update the new solution  

1 k  x
, which is the center of the newly updated trust region, and the radius   1 k   of the trust region for next iteration. The scheme for these updating are as following:

 If   12 01 k r      ,
    1 kk    .  If   1 k r   ,
which implies that the approximate model is not appropriate in the region, thus, we do not accept the solution    x and the radius of the trust region is reduced by a factor [START_REF] Arora | Optimization: Algorithms and Applications[END_REF]

1   . That is,     1 kk   xx ,     1 1 kk      .  If   2 k r
  , which denotes that the approximate model is very appropriate, correspondingly, we do accept the solution    x and enlarge the trust region by a factor 2 1   , i.e.,    

1 k   xx and     1 2 kk      .
Typical values for constants in the algorithm are 1 0.25   , 2 0.75   , 1 0.5   and 2 3.   According to above described procedure, the trust-region moves and updates iteratively until the optimality is found or a fixed number of iterations is reached, which is presented in Algorithm 1.5.

construct an approximate model     k  x in Equation (1.16) for the objective function   f x in current trust-region         k k k d       x x x . 4: find the minimum of     k  x in   k  ,         arg min k k     x xx . 5:
compute the ratio

                      k k k k k r f f              x x x x . 6: if   1 k r   then 7:
reject the step and reduce the trust-region:

        11 1
,

k k k k      xx 8: else if   2 k r   then 9:
accept the step and keep the trust-region:

        11 , k k k      xx 10: else if   2 k r   then 11:
accept the step and enlarge the trust-region : EAs use population-based search mechanism. The general framework of an evolutionary algorithm is presented in Algorithm 1.6.In the initialization of an EA, the first generation (parent population), consisting of one or more individuals, is created, and the fitness values (the objective function value is called fitness in EAs) of its individuals are evaluated. Then, the so-called evolution loop is entered, which consists of the evolution operators recombination (also known as crossover), mutation, evaluation, and selection [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF].

        11 2 , k k k        xx 12: end if 13: 1 kk  14: until     k f   x 1.2.
Recombination creates new individuals, also called offspring, from the parent population.

Mutation adds random changes to the newly created offspring. Subsequently, the fitness of the offspring is evaluated. Based on the fitness, selection identifies a subset of individuals which form the new parent population for the next iteration of the evolution loop. The loop is stopped when the termination criterion is fulfilled.

Other Derivative-free Methods

In addition to the Nelder-Mead method, Trust-region method and EAs that we introduced above, there are numerous other derivative-free optimization methods. As deterministic methods, there are direction search methods including Generalized pattern search method [START_REF] Torczon | On the convergence of pattern search algorithms[END_REF], Mesh adaptive direct search method [START_REF] Audet | Mesh Adaptive Direct Search Algorithms for Constrained Optimization[END_REF], Lipschitzian-based methods, such as DIRECT algorithm [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF] and Brand-and-bound search, and multilevel coordinate search [START_REF] Huyer | Global Optimization by Multilevel Coordinate Search[END_REF], etc. For methods that mimic natural processes or some other physical analogies, besides EAs, there are simulated annealing (SA) algorithm, particle swarm optimization (PSO) algorithm, ant colony optimization (ACO) algorithm, harmony search (HS) algorithm, and so forth.

Review of Surrogate-Assisted Evolutionary Optimization

Surrogate-assisted evolutionary algorithms are mainly motivated from reducing computational cost in evolutionary optimization of expensive problems. In surrogateassisted evolutionary algorithms, the so-called surrogate models or metamodels are constructed to simulate the behavior of the original (true) fitness function. Because surrogate models are much cheaper to evaluate, the computational costs are reduced with the assist of surrogate models. In recent years, surrogate-assisted evolution methods based on fitness approximation are preferable and popular in real-world applications, especially in expensive optimization problems. Mainly based on the review papers [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF][START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF][START_REF] Jin | Surrogate-assisted evolutionary computation: Recent advances and future challenges[END_REF] and book [START_REF] Tenne | Computational Intelligence in Expensive Optimization Problems[END_REF], this section briefly reviews surrogate-assisted evolutionary optimization based on fitness approximation.

Three main aspects of surrogate-assisted evolutionary optimization are: types of fitness approximation methods, the working styles and the management schemes of the fitness approximation [START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF]. The following of this section expands these three aspects successively.

Fitness Approximation Methods

The popular and most commonly used fitness approximation methods are using socalled surrogate models or metamodels that are constructed based on a set of evaluated points from the evaluation history. Based on machine learning and statistical learning techniques, so far, several models have been used for fitness approximation. The most popular ones including polynomial regression, Kriging model, radial basis functions, neural networks and support vector machines. A comprehensive description of surrogate modeling techniques are presented in Chapter 3.

Working Styles of Fitness Approximation

The working style of fitness approximation denotes the mechanism of incorporating the fitness approximation models (surrogate models) into evolutionary algorithms (EAs).

Surrogate models can be embedded in almost every operations of evolutionary algorithms, such as initialization, mutation, recombination and fitness evaluations [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF]. According to [START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF], the incorporation mechanisms of surrogate in EAs can be divided into direct and indirect fitness replacement methods, i.e., direct and indirect styles, as shown in Figure 1.2.

The direct fitness replacement method is to use the approximate fitness to directly replace the original (exact) fitness during the evolutionary optimization. Individuals mostly have the approximate fitness during the optimization. Only a few individuals are evaluated by original fitness function for control purpose. The indirect fitness replacement method is to use the approximate fitness only for some but not all processes in the EAs, such as population initialization and EA operators. The original fitness is kept for each individual and the approximate fitness is not used to directly replace the original fitness.

Direct Style

Direct style, i.e., direct fitness replacement method, is a straightforward strategy to use fitness approximation models. Individuals are evaluated by surrogate models and then the estimated fitness is assigned to each individual. In other words, the fitness approximation model (surrogate model) undertakes the role of the original fitness function, and thus the original fitness is replaced by approximate fitness. However, the obvious drawback is that the evolutionary algorithm may be misled to false optimum due to the inaccuracy of the surrogates of the original fitness function [START_REF] Jin | On evolutionary optimization with approximate fitness functions[END_REF]. A false optimum is an optimum of the surrogate model, which is not an optimum of the original fitness function [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF]. Therefore, in most cases, surrogates should be used together with the original fitness function in order to prevent the optimization process from being misled by false optima introduced by surrogates.

This is termed as model management [START_REF] Dennis | Managing Approximation Models in Optimization[END_REF][START_REF] Yang | Managing approximation models in multiobjective optimization[END_REF] or evolution control [START_REF] Jin | On evolutionary optimization with approximate fitness functions[END_REF][START_REF] Grä Ning | Individual-based Management of Meta-models for Evolutionary Optimization with Application to Three-Dimensional Blade Optimization[END_REF], which is of significant importance in direct fitness replacement methods. The details of model management will be described below in Section 1.3.3.

Indirect Style

In indirect style, i.e., indirect fitness replacement method, the exact fitness of each individual is computed during EA process and the approximate fitness is used in other ways.

For example, the approximate fitness can be used for population pre-selection. In [START_REF] Ulmer | Optimization by Gaussian Processes assisted Evolution Strategies[END_REF], the Gaussian process model was used to pre-select the most promising solutions, which were then actually evaluated by the original fitness function. In this method, based on a standard According to [START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF][START_REF] Rasheed | Comparison of methods for developing dynamic reduced models for design optimization[END_REF][START_REF] Rasheed | Informed operators: Speeding up Genetic-algorithm-based design Optimization using Reduced Models[END_REF], indirect fitness replacement method uses approximate fitness in mutation and recombination (crossover) operators through a technique known as Informed Operators. In this approach, the approximate models are used to evaluated the fitness of candidates only during the mutation and/or recombination process. After the mutation and/or mutation process, the exact fitness is still computed for each individuals. The advantage is that using the approximate fitness indirectly in the form of informed operator rather than direct fitness replacement is expected to keep the optimization moving toward the true global optima and to reduce the risk of convergence to suboptimal solutions or false optima because each individual in the population still carries its exact fitness [START_REF] Rasheed | Informed operators: Speeding up Genetic-algorithm-based design Optimization using Reduced Models[END_REF]. Some of the informed operators in [START_REF] Rasheed | Comparison of methods for developing dynamic reduced models for design optimization[END_REF][START_REF] Rasheed | Informed operators: Speeding up Genetic-algorithm-based design Optimization using Reduced Models[END_REF] are described as follows:

 Informed initialization: An individual in the initial population is generated by selecting the best individual from a number of uniformly distributed individuals in the design space according to the approximate fitness.

 Informed mutation: The informed mutation is best mutation from several random mutations according to the approximate fitness. Specifically, several random mutations of the base point are firstly generated ; these mutations then are evaluated by the surrogate model; the mutation with the best approximate fitness is returned as the result.

 Informed recombination: Two parents are selected randomly according to the usual selection strategy. These two parents are not changed in the course of the informed recombination (crossover) operation. Several recombinations are conducted by randomly selecting a crossover method, randomly selecting its internal parameters and applying it to the two parents to generate a potential child.

The surrogate is used to evaluate every potential child, and the best child is selected as the outcome.

Model Management

The direct working style of fitness approximation uses surrogate model for fitness evaluations and may reduce the number of fitness evaluations significantly [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF]. However, the application of surrogate models to evolution computation is not as straightforward as one may expect. Apparently, there are two principles in direct fitness replacement strategy. First, it should be ensured that the evolutionary algorithm converges to the global optimum or a near-optimum of the original fitness function. Second, the computational cost should be reduced as much as possible. In other words, the number of evaluating the expensive original fitness function should be decreased as much as possible.

Due to the lack of data and the high dimensionality of design space, it is very difficult to construct a perfect global approximate model of the original fitness function. The inaccuracy of the approximate model may lead the EA to inferior local optimum or false optimum. To tackle this problem, two main aspects are considered. Firstly, the approximate model should be used together with the original fitness function. This is known as evolution control in evolutionary computation [START_REF] Jin | On evolutionary optimization with approximate fitness functions[END_REF]. Secondly, the quality of surrogate model should be improved as much as possible with the given limited data. Several aspects are important to improve the model quality, such as selection of the model, selection of surrogate model training method and selection of error measures [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF]. Therefore, model management, which consists of above two aspects, is essential for direct fitness replacement methods.

Evolution Control

Evolution control means that, in surrogate-assisted evolutionary computation, the original fitness function is used to evaluate some/all individuals in some/all generations [START_REF] Jin | On evolutionary optimization with approximate fitness functions[END_REF].

An individual that is evaluated using original fitness function is called a controlled individual.

Similarly, a generation in which all its individuals are evaluated using the original fitness function is called a controlled generation [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF]. The evolution control in surrogate-assisted evolutionary computation generally can be divided into three main approaches [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF][START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF].

No Evolution Control

When the built surrogate model is assumed of high-fidelity, the original fitness function is not used in evolutionary computation, i.e., no individual or generation is controlled. This is known as no evolution control. In this case, after the surrogate model is constructed, the original fitness function is not at all used in evolutionary computation.

Fixed Evolution Control

Fixed evolution control implies that the frequency of evolution control is fixed. There The main drawback of fixed evolution control methods is that the frequency of evolution control is fixed. This is not very practical because the fidelity of the surrogate model may vary significantly during optimization process. As a matter of fact, a predefined fixed evolution control frequency may cause strong oscillation during optimization due to large model error [START_REF] Ratle | Accelerating the convergence of evolutionary algorithms by fitness landscape approximation[END_REF].

Adaptive Evolution Control

It is intuitive and reasonable to assume that the frequency of evolution control should depend on the fidelity of the surrogate model. This strategy is called adaptive evolution control. In adaptive evolution control, individuals are iteratively controlled to update the surrogate model until the fidelity (or quality) of the surrogate model is acceptable.

Off-line Model Training

Off-line model training denotes the training process before the model is used in evolutionary computation [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF]. Off-line model training constructs surrogate models based on data sampling (design of experiments) or previous optimization history data. In this case, either the surrogate model is of high fidelity or the original fitness function is expensive to evaluate, so that the original fitness is not used after the model training [START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF].

On-line Model Updating

On-line model training (updating) denotes rebuilding or updating the model during the evolutionary process. In surrogate-assisted evolutionary computation, surrogate model may be constructed at an early stage of the EA process. Because the set of sample points, which are used to train the model, is limited and does not cover the whole search space, the surrogate may concentrate on the region spanned by the existing sample points and not cover the rest of the search space well. As the EA continues and new individuals enter into the population, the accuracy of the previously built surrogate model will decrease. By evolution control, new points are added into the training data set. Thus the surrogate model needs to be retrained using the old sample points together with the new sample points. This is the typical on-line model updating technique in surrogate-assisted evolutionary computation.

Introduction

Evolution strategies (ESs) are stochastic, derivative-free optimization methods. It has been proven that ESs are appropriate and successful for continuous black-box optimization, i.e., for optimization scenarios, where no analytical expressions of the objective functions are explicitly given and derivatives are unavailable [START_REF] Kramer | Machine Learning for Evolution Strategies[END_REF]. Additionally, ESs have the capability of accessing the global optimum for multimodal problems [START_REF] Kramer | A Brief Introduction to Continuous Evolutionary Optimization[END_REF]. Therefore, ESs are increasingly popular in solving real-world optimization problems.

ESs are search paradigms inspired by the principles of biological evolution. They belong to the family of evolutionary algorithms (EAs) that address optimization problems by implementing a repeated process of stochastic variations followed by selection. In each generation (iteration), new offspring (candidate solutions) are generated from their parents (candidate solutions that have already been visited), their fitness are evaluated, and then the better offspring are selected to become the parents for the next generation [START_REF] Hansen | Evolution Strategies[END_REF].

ESs most commonly address the problem of continuous black-box optimization. The search space is the continuous domain, d , and solutions in search space are d-dimensional vectors, denoted as

d  x
. The objective or fitness function

 

:, d ff  xx is considered to be minimized. In the field of evolutionary algorithms, the objective function is often called fitness function, correspondingly, the objective function value is called the fitness value. High fitness means low fitness function values in the convention minimizing problem. There is no specific assumptions on fitness function f , except that the fitness function f can be evaluated for each x . This kind of search problem is referred to as black- box optimization. The objective is to generate solutions ( x vector) with small fitness values while using a small number of fitness function evaluations.

Main Principles and Evolutionary Operators

Inspired by the principles of biological evolution, evolution strategies use a repeated process that consist of evolutionary operators which mimics the mechanisms of the Darwinian theory of evolution to address optimization problems. A population consists of the so-called individuals is undergoes evolution cycles during each generation (iteration) of evolutionary search. Each individual consists of a candidate solution or vectors of input

parameters d  x
, an associated fitness value   f

x , and of the so-called endogenous parameters, which are strategy parameters for the mutation operator [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF]. In some cases, the population contains only one individuals. Individuals are also denoted as parents or offspring, depending on the context. Generally, a generation loop can be described as :

1. One or several parents are selected from the population (mating selection) and new offspring are generated by recombination of these parents.

2. The new offspring undergo mutation and become new members of the population.

Environmental selection reduces the population to its original size by choosing better individuals.

The sequence of generations is continued until a termination criterion is met. Typical termination criteria are set as reaching a maximum number of evaluations, reaching a target fitness value, or stagnation of the search process.

Evolution Strategies 29

The main principles within above generational procedure of evolution strategies are firstly presented below. Then the corresponding evolutionary operators, which are implementations of main principles in algorithms, used in evolution strategies are described.

Main Principles

Environmental Selection

In evolution strategies, environmental selection is applied as so-called truncation selection. Specifically, based on the individuals' fitness   f x , only  best individuals from the population survive. In other words, only  individuals that have best fitness are selected and kept in the population. In contrast to roulette wheel selection in genetic algorithms [START_REF] Man | Genetic Algorithms[END_REF], only fitness ranks are used in truncation selection. In evolution strategies, environmental selection is deterministic. Environmental selection, on average, increases the fitness of the population and at the same time reduces diversity as some individuals are discarded. Overage individuals can also removed by environmental selection.

Mating Selection and Recombination

In biological filed, recombination, also known as crossover, mixes the genetic material of parents. Similarly, in evolution strategies, recombination combines information from several parents to generate a new offspring. Firstly, mating selection picks individuals from the population to become new parents. These parents are then used to generates a single new offspring by recombination. There are two common scenarios for mating selection and recombination [START_REF] Hansen | Evolution Strategies[END_REF]:  Fitness-independent mating selection and recombination do not depend on the fitness values of the individuals and can be either deterministic or stochastic. In this scenario, environmental selection is essential to drive the evolution toward better solutions.

 Fitness-based mating selection and recombination, where the recombination operator utilizes the fitness ranking of the parents. Thus, recombination performs in a deterministic way. Environmental selection can be potentially be omitted in this case.

Mutation and Parameter Control

Mutation introduces small, random, and unbiased changes (perturbations) to an individual. This provides the main source of variation (evolutionary changes) of individuals in evolution strategies. Thus, mutation plays an important role in evolution strategies. The average size of these changes introduced by mutation, depends on endogenous parameters that change over time, for example, the step-size  and covariance matrix C . These parameters are also called control parameters, or endogenous strategy parameters. In contrast, exogenous strategy parameters are fixed once and for all, for example, the parent number  .

Parameter control, i.e., how to control the endogenous parameters, is not always directly inspired by biological evolution, but it is an indispensable and central feature of evolution strategies.

Evolutionary Operators in Evolution Strategies

In order to realize algorithms of evolution strategies, principles described in above subsection need to be specified in concrete evolutionary operators so that to form evolutionary search schemes. The selection, recombination and mutation operators used in evolution strategies are explained below.

Selection Operators

Selection operator, which relates directly to the environmental selection concept of survival of the fittest, gives the evolutionary search a direction. In each generation, the selection operator emphasizes better solutions by selecting a subset of individuals based on their fitness to form the new parent population used in the next generation. Usually, we use  to denotes the number of individuals in parent population (also called the size of parent population), and  to denotes the number of offspring generated in each generation. In Neglecting superior solutions, in comma selection, may sound irrational. However, good solutions can be only local optimum. By using comma selection, the evolutionary process could escape from the local optimum and reach a better optimum [START_REF] Kramer | Machine Learning for Evolution Strategies[END_REF]. In contrast, when using plus selection, the search process may fail to leave the local optimum without the ability to neglect [START_REF] Kramer | A Brief Introduction to Continuous Evolutionary Optimization[END_REF]. Thus, comma selection is advantageous in the case of multimodal problems.

Recombination Operators

In In evolution strategies, the intermediate or weighted recombination with   is often used. Thus, the result of selection and recombination is often deterministic. This means that eventually all offspring are generated by mutation from the same single solution vector (the parent centroid, generated by recombination).

Mutation Operators

The mutation operator, which introduces variations by adding perturbation to the result of recombination, provides the main source of variation (evolutionary changes) of offspring in evolution strategies. According to [START_REF] Vanneschi | Handbook of Natural Computing[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF], a mutation operator is supposed to fulfill three rules, namely:

 Reachability: from an arbitrary solution, any point in the search space can be reachable with probability strictly larger than zero by means of a finite number of applications of the mutation operator.

 Unbiasedness: the mutation operator should be unbiased, unless knowledge about the problem has been gathered. This can be achieved by using the maximum entropy distribution that obeys knowledge about the problem as constraints.

 Control: the mutation operator should have parameters that affect the mutation strength (shape of the distribution), such that the extent of variation can be controlled. As known from theory, when approaching the optimal solution, the strength of mutation must be weakened steadily.

In continuous search space d , these principles are fulfilled by the famous Gaussian mutation operator [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF], which uses the multivariate normally distributed random vector as the

 

, d mC is sometimes used like a vector which is actually sampled according the distribution given. Since the covariance matrix C is symmetric and positive definite, the following eigendecomposition exists: where B is an orthogonal matrix (i.e., TT  B B BB I where I is identity matrix), the columns of which are the eigenvectors of C , D is a diagonal matrix with square roots of eigenvalues of C as diagonal elements, and

1 2 dd   C
is the symmetric square root matrix of C such that

  1 1 1 1 2 2 2 2 T  C C C C C .
With above eigendecomposition and considering the properties of multivariate normal random vector [START_REF] Rencher | Methods of multivariate analysis[END_REF], the multivariate normal distribution   , d mC can be written as

      1 2 ,, ,, dd d 
 m C m 0 C m C 0 I (2.6) where 1 2
T  C BDB is given in Equation (2.5).

In evolution strategies, the Gaussian mutation operator generates a mutated individual

 x from an individual d  x
by adding a multivariate normally distributed vector on x . This can be expressed generally as

    1 2 , , . dd      x x 0 C x C 0 I (2.7)
Based on multivariate normal distribution, the three different mutation operators that are commonly used in evolution strategies are:

(a) Isotropic or Spherical Gaussian mutation generates the mutation  x from x by using  BI and   DI for matrices B and D in Equation (2.7), where  is called step-size, that is,

  ,. d    x x 0 I (2.8)
This corresponds with spheres with individual radii defined by  , as indicated in (c) Correlated Gaussian mutation denotes the situation that matrix B is not just an identity matrix and

    1 diag diag , , d   D  is a diagonal matrix, i.e.,           2 diag , , diag , d d d        x x B 0 I x B 0 x 0 C  (2.10)
This rotates the hyperellipsoids with respect to the coordinate axes, as shown in the right part of Figure 2.1 (c).

The choice of the mutation operator from above described three cases has a direct impact on the complexity of the endogenous parameters that control the multivariate normal

Parameter Control and Algorithms

Controlling the parameters of mutation operator is essential to the design of evolution strategies. Consider an evolution strategy using the isotropic mutation operator, where the step-size  is a scaling factor for the random vector perturbation. The step-size  controls to a large extent the convergence speed. In situations where larger step-sizes lead to larger expected improvements, a step-size control technique should aim at increasing the step-size (and decreasing it in the opposite scenario) [START_REF] Hansen | Evolution Strategies[END_REF]. Generally, the goal of parameter control is to drive the endogenous strategy parameters close to their optimal values. The research on ESs mainly focuses on question of adaptation of parameters of mutation operators. Different parameter control strategies bring about different algorithms for evolution strategies. In the following, three specific ESs are outlined, each of them representing an important achievement in parameter control.

The 1/5th Success Rule

The 1/5th success rule is a basic step-size control strategy for evolution strategies. It is based on an important discovery made by Rechenberg [36] in early research of evolution strategies. The (1+1)-ES using 1/5th success rule is the first evolution strategy and the basis of the evolutionary algorithm in the field of ES. This algorithm is implemented in Algorithm 2.1.

As the first and simplest evolution strategy, (1+1)-ES only use mutation and selection operators to guide the search. We denotes a generation (iteration) of the search by the generation counter g  ( 0,1, 2, g 

). In generation g , a single offspring In (1+1)-ES, the adaptation of the step-size, which is the only one strategy parameter of (1+1)-ES, is of crucial importance for reliable results and efficiency of the algorithm. The step-size significantly affects the convergence speed. The step-size has to be adapted in order to speed up the optimization process. Usually, in situations where larger step-size lead to larger expected improvements, a step-size control technique should aim at increasing the step-size, and decrease it in the opposite situations. The step-size  in (1+1)-ES is adapted according to the 1/5th success rule.

The idea of the 1/5th success rule is to increase the mutation rate (step-size in isotropic mutation), if the success probability (i.e., the empirical probability that  

    1 g ff  xm
), which refers to the ratio between successful mutations and all mutations, is larger, and to decrease it, if the success probability is smaller. Specifically, if about 1/5 of all mutations are successful, the step-size is optimal and no adaptation is required. The step-size needs to be reduced, if the success rate is smaller than 1/5, and needs to be increased, if the success rate is larger than 1/5. The new step-size   1 g   is updated according to Line 10 in Algorithm 2.1. In the algorithm, the empirical probability that  

    1 g ff  xm is expressed by         1 g ff     xm
, where    denotes the expected value function, and    is the

indicator function such that         1 1 g ff  xm only if       1 g ff  xm and         1 0 g ff  xm
otherwise. The success probability information is obtained from previous iterations.

Cumulative Step-Size Adaptation (CSA)

The Cumulative Step-Size Adaptation Evolution Strategy (CSA-ES) developed by Ostermeier et al. [START_REF] Ostermeier | Step-size adaptation based on non-local use of selection information[END_REF] adapts the step-size  by using a so-called evolution path that records information accumulated from the preceding generations instead of using information from the current generation only. The evolution path records the sum of consecutive successful mutation steps to make a decision about possible corrections of the step-size. This information can improve the adaptation and search procedure remarkably [START_REF] Hansen | Evolution Strategies[END_REF]. When successful steps are towards the same direction, the evolution path which records the sum of mutation steps will be relatively long in this direction. Because the same distance in this direction can be covered with larger steps, the step-size should be increased [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF]. On the contrary, if the orientations of successful steps are opposite to each other, they sum up making the evolution path relatively short. In this case, the used step-size is too large and should be decreased. In this algorithm, comma selection, isotropic mutation and intermediate recombination are used. The step-size  is adapted by using the evolution path  p . In the evolution loop, at generation g ,  offspring   1 ,,  xx are generated by adding a isotropic Gaussian mutation on the solution vector   g m (also known as the mean of mutation distribution), which is described in Lines 4~6. The fitness of these  individuals then are evaluated (Line 7). The selection operator chooses  best individuals to form the new parent population according to their fitness values. Then, the intermediate recombination of the  parents generates a single new solution vector  

1 g  m
, which is the new mean vector of the mutation distribution (Line 9).

The evolution path  p is updated (Line 10) using local information about a successful mutation step 1

: 1 i i      z , a
  , 4 3ln , 2 dd              ,           2 5 , 1 2 max 0, 1 1 1 . c d d d c                 2: initialize       0 0 0 , 0, , 0 d g       m p 0 3: repeat 4: for 1, , k   do 5:   , kd  z 0 I //i.i.d. for each k z 6:     gg kk   x m z //mutation 7:   kk ff  x 8:
end for 9:

  1 1 : 1 g i i        mx //

Mutation or Sampling

In  

, -CMA-ES w    , the  offspring are generated by correlated mutation from the same single solution vector, i.e.,

                1 2 , , for 1, ,
,

g g g g g g k d k k         x m 0 C m C z (2.11)
where: 0 g  is the counter of generation.

 is the number of offspring. Its default value is

  4 3ln d     ( d is the dimension of the problem). d k 
x is the k-th offspring (individual, search point, solution).

  g d  m
is the mean of mutation distribution (a solution vector), which is formed by weighted recombination of selected  best individuals after initialization. Its initial value   0 m is given by the user.  is given according to the bounds of search space.

  g dd   C
is a symmetric positive definite covariance matrix at generation g . Its

default initial value is   0  CI .     1 2 t dd   C is the symmetric square root matrix of   g C such that                   1 1 1 1 2 2 2 2 T g g g g g     C C C C C
, which can be obtained according to Equation (2.5)

d k  z are random vectors (mutation vector) sampled from   , d 0I , i.e.   , kd  z 0 I .
From Equation (2.11), the  offspring also can be viewed as sampling from a multivariate normal distribution with mean vector  

g m and covariance matrix       2 gg  C , i.e.,           2 , g g g kd  x m C .
To define the complete iteration procedure, the remaining question is how to calculate  

1 g  m ,   1 g  C and   1 g   for the next generation 1 g  .

Selection and Recombination

After generating  offspring, their fitness function values are evaluated. According to their fitness,  (   ) best individuals from  offspring are selected to become the parent population and undergo weighted recombination to generate a single new solution vector

  1 g  m . Thus, the new solution   1 g  m
is computed as a weighted average of the selected points:

          1 1 2 :: 11 u g g g g i i i i ii ww          m x m C z (2.12) 12 1 1, 0 i i w w w w          (2.13)
where:  is the number of selected individuals. Its default value is

        1 ln 2 1 2 ln ln 2 1 2 ln i j w i j                  , 1, , i   . : i  x signifies the i-th best individual out of 12 , , ,  x x x from Equation (2.

11). The index :

i  denotes the index of the -th i ranked individual and

      1: 2: : f f f        x x x
where f is the objective function.

: i  z is the mutation vector associated with :

i  x .

Adaptation of Step-size and Covariance Matrix

The CMA-ES uses accumulative step-size adaptation (CSA) to adapt the step-size  and the covariance matrix adaptation (CMA) to update the covariance matrix C . These two techniques (CSA and CMA) use so-called evolution paths for accumulating strategy parameter information across several generations [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF]. Two evolution paths, The evolution path  p and the step-size  are recursively computed [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF][START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF] according to the equations below: (2.15) where:

                        1 1 1 2 : 1 12 12 gg g g t w t g w i i i c c c c c c w                             mm p p C pz (2.14)           1 1 exp 1 , g gg d c d                p 0I
  g d   p
is the evolution path for global step-size at generation g , its default initial

value is   0   p0 .   0,1 c  
is the time constant for the adaptation of the step-size, its default value [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF] is

    25 ww cd       .
w  is the so-called variance effective selection mass, which is defined by:

  2 2 11
.

w i i ii ww      If equal weights 1 i w   are used, w  is equal to  .     1 2 g  C is the inverse of     1 2 g C which is the symmetric square root matrix of   g C .
d  is the damping factor, its default value is

1 1 2 max 0, 1 1 w dc d             .     , d 0I is the expectation of the Euclidean norm of a   , d 0I distributed random vector, which can be approximated by     2 11 1 , 2 1 22 4 21 d dd d dd                      
The evolution path c p and covariance matrix C are updated [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF][START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF] according to the equations below:

                        1 1 c c c c c 1 2 c c c c : 1 12 12 gg gg w t gg w i i i c h c c c h c c w                    mm pp p C z (2.16)                           T 1 1 1 1 1 c c T 11 22 :: 1 1 g g g g g u u gg i i i i c c c h cw                             C C p p C C z C z
(2.17)

where:

  c g d  p
is the evolution path for covariance matrix at generation g , its default initial

value is   0 c  p0 .   c 0,1 c 
is the constant for covariance matrix adaptation. c 1 c is the backward time horizon of the evolution path c p . The default value of c c is:

    4 4 2 c w w c d d d     
h  is a Heaviside function defined by:

          1 21 2 1 if 1.4 , 1 11 0 otherwise g d g d h c                    p 0I
The purpose of function h  is to avoid an update of c p to take information of the current generation into account, when

  1 g   p
is too large [START_REF] Hansen | The CMA Evolution Strategy: A Comparing Review[END_REF]. This prevents a too fast increase of axes of C when the step-size is far too small.

     

cc 12 h h c c     
is of minor relevance. In the unusual case of 0 h   , it substitutes for the second summand from Equation (2.16) in (2.17) [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF].

In summary, for CMA-ES iteration, the new offspring population of search points is generated by Equation (2.11), the selection and recombination is defined by Equation (2.12), the adaptation of step-size is given by Equations (2.14) and (2.15), and the covariance matrix is updated with Equations (2.16) and (2.17). The iteration defined above is repeatedly executed until the termination criterion fulfilled. Putting it all together, the pseudocode of  

, w    -CMA-ES is given in Algorithm 2.3. Algorithm 2.3 The   , -CMA-ES w    1: given:   , 4 3ln , 2 dd             ,         1 ln 2 1 2 ln for 1, , ln 2 1 2 ln i j i wi j            ,   2 2 11 w i i ii ww      , 21 , 1 2 max 0, 1 51 ww w c d c dd                       11 22 4 2 1 2 , , min 1 , 2 . 42 1.3 2 w w w cu w ww d c c c c dd dd                      2: initialize           0 0 0 0 0 c , 0, , , , 0 d g 
       m p 0 p 0 C I 3: repeat 4: for 1, , k   do 5:   , kd  z 0 I //i.i.d. for each k z 6:         12 g g g kk   x m C z //mutation 7:   kk ff  x 8:
end for 9:

          1 1 2 : 1 g g g g ii i w        m m C z 10:         1 : 1 12 gg w i i i c c c w                p p z 11:           1 1 exp 1 , g gg d c d                 p 0I 12:             1 1 2 c c c c c : 1 12 g g g w i i i c h c c w           p p C z 13:                           T 11 T 1 1 1 22 1 1 c c : : 1 1 u g g g g g g g u i i i i c c c h c w                     C C p p C C z C z
14:

1 gg  15: until termination criterion is fulfilled

Surrogate Modeling of Computer

Experiments

This chapter reviews the surrogate modeling method, which is an approach to model the (input-output) behavior of the simulation model. Computer simulation tools, such as finite element analysis (FEA) and computational fluid dynamics (CFD), are more and more widely applied in engineering problems. Many engineering design problems require a number of simulations to evaluate design objective. However, for many read world problems, a single simulation can take many minutes, hours, or even days to complete. Consequently, computation-intensive tasks, such as design optimization, sensitivity analysis and reliability analysis, become impractical or impossible because they require doing hundreds, thousands or even millions of simulations. One way of reducing this computational burden is by constructing approximation models, known as surrogate models, or metamodels, that approximate the input-output relationship of the simulation model as closely as possible while being computationally cheap to evaluate [START_REF] Fang | Design And Modeling for Computer Experiments[END_REF]. Then, many tasks can be implemented by using the built surrogate model. Computational cost is thus apparently reduced owing to the aid of cheap-evaluated surrogate model. Main issues about surrogate modeling are presented in the following of this chapter. In

Surrogate Modeling Process

Surrogate modeling is concerned with the construction of mathematical models to describe the relationships between specific inputs and outputs exhibited by the simulation model, based on a set of limited data acquired by running the simulation model with intelligently chosen input variables. Since simulation models usually generate multiple (a set of) responses (outputs), a surrogate model can be constructed for each of the outputs of interest. Without loss of generality, only one response of interest is considered in this chapter.

Let       T 11 , , , , , , d dd y f f x x x x S      xx (3.1)
denotes a simulation model which models a physical system, where

  T 1 ,, d xx  x is a
vector of d design variables or input variables (input), y is a response (output), S is referred to as the design space or input variable space, and the function f signifies the computer simulation model, such as FEA and CFD model of a problem, that maps the d-dimensional inputs into a scalar output. Model (3.1) generally can be regarded as a solution of a set of equations, including linear, nonlinear, ordinary and/or partial differential equations, and it is often impossible to obtain an analytic solution for the equations.

Then, the general surrogate modeling problem can be stated as follows: "Given a set of samples (observations)

    T 1 ,, n    X x
x in the design space and the corresponding response values

            T T

Design of Experiments

Design of experiments (DOE, or experimental design) refers to the techniques that are used for guiding the choice of experiments to be performed in an efficient way. DOE techniques were originally developed to study the behavior of systems through physical experiments [START_REF] Box | Empirical Model-Building and Response Surfaces[END_REF]. The basic objective of DOE is to determine multiple combinations of the controlled parameters (or conditions) at which the experiments will be conducted. Each combination of controlled parameters, in mathematical terms, can be considered a sample.

Although DOE originally referred to physical experiments, in modern times, it would includes both physical and computer experiments. This section only focus on the design of computer experiments.

In computer experiments, DOE is a strategy for allocating samples (points) in the design space that aims at maximizing the amount of information acquired [START_REF] Koziel | Computational Optimization[END_REF]. Computer simulations (experiments) are performed at these points to create the training data set that is subsequently used to construct the surrogate model. Obviously, there is a trade-off between the number of samples and the amount of information that can be extracted from these samples. Thus, a good experimental design should minimize the number of samples needed to acquire enough information.

In following part of this section, the factorial designs, which are usually referred to as classical DOE techniques, and several space-filling or modern DOE techniques are presented.

Before giving descriptions of these DOE techniques, some terminology in design of computer experiments are introduced. The controllable variables that are of interest in computer experiments are usually called design variables or input variables or simply called inputs. Design space or input variable space is the space where the input variables take values.

A point in design space or an input setting, which is a combination of input variables, is called a sample (sample point). A run is the implementation of an experiment (computer simulation) with the input variables given in a sample. Output(s) or response(s) is the result of a run based on the purpose of the experiment.

Factorial Designs

Factorial designs, which are classical DOE techniques [START_REF] Giunta | Overview of modern design of experiments methods for computational simulations[END_REF], are straightforward techniques for experimental design. In factorial designs, the range of each design variable is divided into different levels between the upper and lower limits of the design space. Factorial designs allocate samples at combinations of different levels of all the design variables. In other words, a factorial design is a set of level-combinations of the factors. A factorial design is called symmetric if all factors have the same number of levels ; otherwise, it is called asymmetric.

When samples are located at all the combinations of the different levels of all the design variables, the design is called full factorial design. Obviously, the number of sample points in a full factorial design should be 1 d j j nq    , where j q is the number of levels of the design variable j . When all the design variables have the same number of levels q , the number of samples are d nq  . Figure 3.2 illustrates a full factorial design of three-dimensional design space with three level for all variables. The number of sample points of full factorial design increases exponentially with the number of design variables. Thus, in high-dimensional problems, the full factorial design approach may be cost/time prohibitive. Therefore, we consider implementing a subset of all the level-combinations that have a good representation of the complete combinations. In situations that the running of an experiment is expensive and the number of design variables is large, fractional factorial design (FFD) can be used. Fractional factorial design (FFD) uses a fraction of the full factorial design sample points, i.e., a FFD is a subset of all level-combinations of the design variables. Consequently, how to choose a good subset is the most important issue in FFD. A carefully selected combination known as the orthogonal array is recommended in the literature and has been widely used in practice. In practice, other alternative factorial designs, such as central composite design, star design, and Box-Behnken design, can be used.

Latin Hypercube Designs

This subsection discusses some modern DOE techniques [START_REF] Giunta | Overview of modern design of experiments methods for computational simulations[END_REF] for planning computer experiments, which tend to allocate samples throughout the design space as uniformly as possible. Such designs are broadly referred to as space-filling designs.

Because of the deterministic feature of computer experiments, i.e., samples with the same input setting will yield identical outputs, the basic principles of experimental design for controlling noise and bias, replication, blocking, and randomization, are irrelevant to computer experiments [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF]. The exact functional form of the relationship between inputs and the response is unknown and often very complicated, although the response can be computed at any given inputs. Various surrogate models can be built using different techniques. However, before data are collected, quite often little priori or background knowledge is available about which model would be appropriate, and designs for computer experiments should facilitate diverse modeling methods [START_REF] Dean | Handbook of Design and Analysis of Experiments[END_REF]. For this purpose, a spacefilling design is the best choice. The design region in which to make prediction may be unspecified at the data collection stage. Therefore, it is appropriate to use designs that represent all portions of the design space. When the primary goal of experiments is to make prediction at unsampled points, space-filling designs allow us to build a predictor with better average accuracy.

A variety of space-filling designs are available. One of the most popular category of space-filling designs techniques is the Latin hypercube design, which is based on Latin hypercube sampling. Basic Latin hypercube design and its variants are presented in the following. In order to allocate n sample points in d-dimensional design space with LHS, the range of each variable is firstly divided into n on-overlapping intervals with equal probability, which yields a total number of d n bins in the design space. The samples are randomly selected in the design space so that (i) each sample is randomly placed inside a bin, and (ii) for all one-dimensional projections of the n samples and bins, there is exactly one sample in each bin [START_REF] Koziel | Computational Optimization[END_REF]. Figure 3.3 shows a LHS realization of 10 samples for two design variables. Mathematically, the algorithm for generating LHS samples in design space  

Basic Latin Hypercube Designs

0,1 d is       , 1 , 1 ii i j j j U x i n j d n        (3.2)
where d is the number of design variables, n is the number of samples,    

1 ,, n jj  are independent uniform random permutations of the integers 1 to n , and   i j U are independent random numbers from  0,1 . The subscript j denotes the dimension index and the superscript   i denotes sample number. Then the sampling plan

    T 1 , , , n    X x
x where

      T 1 ,, i i i d xx    x , is a LHD denoted by   , LHD n d . If each ij U in Equation (3.
2) is taken to be 0.5, the corresponding LHS is called midpoint Latin hypercube sampling or centered Latin hypercube sampling. For each variable, Latin hypercube designs have exactly one point in each of the n intervals. This property is referred to as one-dimensional uniformity.

One drawback of Latin hypercube sampling is that there is more than one possible arrangement of bins and samples that meets the LHS criteria. So, a randomly generated Latin hypercube design does not necessarily perform well with respect to criteria of uniformity.

For example, samples allocated along the design space diagonal satisfy conditions of LHS, but they are non-uniformly distributed. Fortunately, there are extensions to the basic LHS technique that can improve the uniformity of Latin hypercube designs. In extensions of basic LHS, the Latin hypercube designs based on distance are of desirable properties.

Latin Hypercube Designs based on Measures of Distance

To construct space-filling Latin hypercube designs, one natural approach is to make use of distance criteria. One of the most widely-used measures to evaluate the uniformity ('space-fillingness') of a sampling plan is the maximin metric introduced by Johnson et al. (1990). The criterion based on this may be defined as follows.

Let   1 i x and   2 i
x be two design points in the design space

  0,1 d D  . For 0 p  , the distance between   1 i x and   2 i
x is defined as

          1 2 1 2 1 1 , . p d p i i i i jj j d x x       xx (3.3)
When 1 p  and 2 p  , the measure in Equation (3.3) becomes the rectangular and Euclidean distance, respectively. Let 1 ,, m dd be the list of unique values of distances between all possible pairs of points in a design X , sorted in ascending order. Further, let 1 ,, m JJ be defined such that j J is the number of pairs of points in X separated by the distance j d .

The simple maximin distance criterion seeks a design X in the design space S that maximizes the smallest inter-point distance 1 d among all available designs. This criterion attempts to place the design points such that no two points are too close to each other.

Clearly, this definition could be applied to any set of sampling plans, but, since we would like to keep the appealing stratification properties of Latin hypercube, we restrict our scope to that class of designs. Nonetheless, even across this narrower domain, Definition 1.1 might still yield several maximin designs. Therefore, we shall use the more complete 'tiebreaker' definition of Morris and Mitchell (1995). Thus, we call design X a maximin design among all available designs if it maximizes 1 d , among designs for which this is true, minimizes 1 J , among designs for which this is true, maximizes 2 d , among designs for which this is true, minimizes 2 J ,…, minimizes m J .

An extended definition of a maximin design was given by Morris and Mitchell (1995).

In order to rank competing designs, based on above definition, they introduced a computationally efficient scalar-value criterion 1 1

q m i q q i i J d       (3.4)
where q is a positive integer. The smaller the value of q  , the better the space-filling properties of X will be. This scalar value distills the cumbersome definition of the maximin criterion into a rather neat and compact form, but it raises the question of how to choose the value of q in Equation (3.4). Values of q are chosen depending on the size of the design searched for, ranging from 5 for small designs to 20 for moderate-sized designs to 50 for large designs [START_REF] Dean | Handbook of Design and Analysis of Experiments[END_REF]. The idea of above described maximin design sounds simple and desirable.

However, generating maximum Latin hypercube designs is a challenging task particularly for large designs.

Other Space-filling Designs

Besides Latin hypercube designs, other space-filling design techniques including orthogonal array, quasi-Monte Carlo sampling, and uniform design can be found in literature.

A good summary of space-filling designs for computer experiments can be found in [START_REF] Giunta | Overview of modern design of experiments methods for computational simulations[END_REF] and Chapter 5 of [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Sampling techniques can be improved by minimizing a specific nonuniformity measure with optimization techniques. The above described maximin Latin hypercube design is a typical instance for this.

Surrogate Models

Once the training data set has been collected from an experiment, the next step is to construct a surrogate model which describes the relationship between the inputs and output(s). Generally, the n samples in d-dimensional design space from DOE can be represented by a nd  matrix X and their corresponding responses are expressed by a 1 n  vector y as

                        1 1 1 1 1 1 1 , d n n n n n d f x x y x x y f                                         x x Xy x x (3.5)
where   , 1, ,

i in  x
denotes the i-th sample, which is a combination of d design variables, and   , 1, , , 1, ,

i j x i n j d 
denotes the j-th component of the i-th sample, and function f signifies the computer simulation model that can calculate the response with any given inputs,

  , 1, , i y i n 
signifies the corresponding response with inputs

  i x . With the training data set                 11 , , , , , nn yy  X y x x
, our goal is to construct a cheap-to-evaluate surrogate model f that emulates the response of expensive computer model f and thus we can use the surrogate model to predict the output of the simulation model for any untried site

x . In this section, we describe several popular surrogate modeling techniques.

Polynomial Regression Model (PR)

Polynomial regression (PR) models have been widely used by researchers for modeling computer experiments. PR model uses a polynomial form model to approximate the inputoutput relationship. The number of terms included in a PR model depends on the desired accuracy of the surrogate and the exact functional form of the relationship between inputs and the response [START_REF] Huang | Uncertainty analysis of deep drawing using surrogate model based probabilistic method[END_REF][START_REF] Messac | Optimization in Practice with MATLAB®: For Engineering Students and Professionals[END_REF]. The general PR model with p basis functions can be expressed as

      T 1 ˆp jj j fb     x x b x  (3.6) where , 1, , j jp   are regression coefficients,  , 1, , j b j p  x are basis functions of the PR model, T 1 ,, p     
is the vector of regression coefficients and

      T 1 ,, p bb    b x x x
is the vector of basis functions. Given the training data set where ŷ is the vector of responses given by PR model and B is the so-called design matrix or Vandermonde matrix, which are respectively defined by Equation (3.8) and (3.9).

                11 , , ,
            T T 11 ˆˆˆ, , , , nn f f y y      y x x (3.8)                         1 1 1 1 1 p n n n p bb bb                          b x x x B b x x x (3.9)
A commonly used approach to estimate the regression coefficients  is the least squares method (LSM). The least squares approach is to minimize the sum of squared residuals at the points in the training data set [START_REF] Fang | Design And Modeling for Computer Experiments[END_REF], i.e., to minimize

            2 2 2 1 1 1 ˆ. p nn i i i i jj i i j S y y y b                 x y B (3.10)
If np  and matrix B is full rank, the least squares estimator of regression coefficients  is given by

  1 TT ˆ.   B B B y  (3.11) Let   T 1 ,, d xx  x
be the d-dimensional design variables, the p polynomial basis functions of first-order (linear) and second-order (quadratic) PR models are:

Linear, 1 pd  :       1 2 1 1 1, , , dd b b x b x     x x x (3.12) Quadratic,    1 2 12 p d d    :                     1 2 1 1 2 2 1 3 1 2 2 1 1 2 2 2 2 3 2 3 2 1 2 2 1 ,, , , , , , , dd k
d d d k d d d pd b b x b x b x b x x b x x b x b x x b x x bx                 x xx x x x x x x x (3.13)
The number of polynomial basis functions dramatically increases with the number of input variables and the order of the polynomial. Thus, lower-order polynomials such as first-

Kriging Model (KG)

Kriging is a Gaussian process [START_REF] Seeger | Gaussian processes for machine learning[END_REF] based modeling method to interpolate deterministic noise-free data [START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF] and has proven to be useful in a wide variety of fields [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF][START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF]. Kriging model is also referred to as "Design and Analysis of Computer Experiments" (DACE) model [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF]. The distinctive feature of Kriging model is that it provides not only a predicted response at an unsampled point but also an estimate of the prediction variance. This variance gives indication of the uncertainty in the Kriging model, which results from the construction of the covariance function. The covariance function is based on the idea that when input points are near one another, the correlation between their corresponding outputs will be high. As a result, the uncertainty associated with the model's predictions will be small for input points which are near the training points, and will increase as one moves further from the training points.

The Kriging model consists of two components: a regression function (also known as global trend function)   g x which is constructed based on the training data set, and a

Gaussian process   Z x which accounts for the local deviation of the data from the trend function. Thus, the general form of Kriging surrogate model is given as

      f g Z  x x x (3.14)
where   g x is a regression function and   Z x is assumed to be a realization of Gaussian process with zero mean, variance 2  and correlation matrix  . The idea is that regression function captures the largest variance in the data (the general trend) and that the Gaussian process interpolates the residuals. The covariance matrix of

  Z x is defined as       2 Cov , , Z Z R      x x x x (3.16)
where

    Cov , ZZ    xx is the covariance of  
                11 , , , , , nn yy  X y x x
, where

    T 1 ,, n    X x
x are n samples in d-dimensional design space, and

    T 1 ,, n yy    y are
the corresponding responses. The prediction mean and prediction variance of Kriging [START_REF] Sacks | Design and Experiments of Computer Experiments[END_REF] can be respectively given by

        T T 1 T f      x b x r x R y F β β (3.18)           T T 2 2 T T ˆT 0 1 f s               bx F x b x r x rx FR (3.19) where       T 1 ,, p bb    b x x
x is the regression basis function vector of the predicting

point x , T 1 ,, p     
is the regression coefficients,   rx is the vector of correlation functions between the untried point and the n sampled (observed) points (defined in Equation (3.20)), R is the nn  matrix of correlation functions for the sampled data (defined by Equation (3.21)), and F is the np  model matrix defined by (3.22). Since the Kriging method assumes that the observed responses are from a Gaussian process, the responses at sampling sites are considered to be correlated random functions with the corresponding likelihood function given by

              T 12 , , , n R R R    r x x x x x x x (3.20)                   (1) (1) (1) (2) (1) ( ) (2) (1) (2) (2) (2) ( ) ( ) (1) ( ) (2) ( ) ( ) , , , , , , , , , n n n
n n n R R R R R R R R R          x x x x x x x x x x x x R x x x x x x (3.21)                 11
        T T 1 T 2 /2 2 1/2 2 1 , , exp 2 2 n L           S y F R y F R β β β θ (3.23)
To simplify above likelihood function, we take the natural logarithm and gives

          T T 1 T 2 2 1 ln ln 2 ln 2 2 2 2 nn L          y F R y F R β β (3.24)
By taking derivatives of Equation (3.24) with respect to β and 2

 respectively, and setting to zero, we obtain maximum likelihood estimates (MLEs) for  and 2   

1 T 1 T 1 ˆ   F R F F R y β (3.25)     T 2 T 1 T 1 ˆn      y F R y F β β ( 3
        2 1 ln ln 2 ln ln 2 2 2 nn LL       R θ (3.27)
The estimator θ of parameter vector θ is obtained by maximize Equation (3.27) under the constraint 0, 1, 2, , l ld   . In other words, parameter vector θ is achieved by solving the following optimization problem:

  Max s.t. 0, 1, 2, , l L lk        θ θ (3.28)
This optimization problem can be carried out using numerical optimization technique.

A global search method such as a genetic algorithm or simulated annealing usually produce good result. After getting θ , the estimators β and 2  can be calculated by Equation

Radial Basis Function Model (RBF)

Radial basis function (RBF) model, which is also called Radial basis function network, can be interpreted as an artificial neural network that uses radial basis functions as activation functions, as shown in Figure 3.4. A RBF network consists of three layers: an input layer, a hidden layer with a non-linear RBF activation function, and a linear output layer [START_REF] Wu | Using Radial Basis Function Networks for Function Approximation and Classification[END_REF].

Usually, the neurons (nodes) in the hidden layer use the RBF in form of Commonly used non-parametric basis functions are [START_REF] Messac | Optimization in Practice with MATLAB®: For Engineering Students and Professionals[END_REF] (writing r  xc ):

      , i i   x x c 1, , im  ,
 Linear:   rr   ,  Cubic:   3 rr   ,
and

 Thin plate spline:     2 ln r r r   .
To improve the generalization properties of RBF model (3.29), parametric basis functions can be adopted. Commonly used types of parametric basis functions include:

 Gaussian:     22 2 r re     ,  Multi-quadric:    
    T 1 ,, n    X x x
and their corresponding responses

    T 1 ,, n yy    y
, the interpolation condition [START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF] can be written in matrix form as  Φwy (3.30) where

                                                      1 1 1 2 1 2 1 2 2 2 12 m m n n n m                             x c x c x c x c x c x c Φ x c x c x c (3.31)
It is easy to see that if mn  and the centers of the radial basis functions coincide with the training points

    , 1, , ii in  cx
, Φ is a regular square matrix and the matrix equation represent highly nonlinear data, which is otherwise often challenging to accomplish using low order PR models [START_REF] Messac | Optimization in Practice with MATLAB®: For Engineering Students and Professionals[END_REF].

Multi-layer Perceptron Networks (MLP)

Multi-layer perceptron (MLP) network, which is a feedforward artificial neural network (ANN) model that maps sets of input data onto a set of appropriate outputs, is a more recent and increasingly popular choice of surrogate model. A MLP network consists of an input layer, several hidden layers, and an output layer. The neuron or node which includes a summation and activation function is the basis structure of ANN.

A simple MLP network (shown in Figure 3.5), that consists of one input layer, one hidden layer with nonlinear activation function, and one output layer with linear activation function, approximates the inputs and outputs as follows:

  00 11 ˆmd j ij i j ji f g w x w          x (3.32)
where

  T 1 ,, d xx  x
is the input, m is the number of nodes (neurons) of the hidden layer, j  is the weight connect between the output and the j-th note in the hidden layer, the function

 

g  is the activation function of the hidden layer, ij w is the weight connection between the i-th component of the input and the j-th note in the hidden layer. The most commonly used activation function is the logistic function (or called Sigmoid function), which has the form:

    1 1 exp g net net      (3.33)
where  is a constant and 0 1

d ii i net w x w   
is the so-called net input signal to the neuron.

A parameters include the number of neurons, the number of layers, the type of activation function, and the optimization algorithm used to train the network. In addition, the training process generally needs to be supervised in order to avoid "over-fitting" [START_REF] Messac | Optimization in Practice with MATLAB®: For Engineering Students and Professionals[END_REF].

Support Vector Regression (SVR)

Support vector machines (SVM) are popular models for solving classification and regression problems in recent years. An important and favorable property of SVM models is that the determination of the model parameters corresponds to a convex optimization problem [START_REF] Montaño | A Review of Techniques for Handling Expensive Functions in Evolutionary Multi-Objective Optimization[END_REF], i.e., there is no local minima during training SVM models, and thus the optimization process does not depend on the problem dimensions and overfitting is seldom an issue [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF][START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF]. In support vector regression (SVR), the goal is to construct a model   f x that has at most an  deviation from the actual targets   i y for all the training data, and at the same time is as flat as possible [START_REF] Smola | A tutorial on support vector regression[END_REF]. In other words, the errors are not considered if they are less than  , but any deviation larger than this will not be accepted. ) is feasible [START_REF] Forrester | Engineering Design via Surrogate Modelling: A Practical Guide[END_REF]. However, sometimes such a solution may not actually exist and it is also likely that better predictions will be obtained if we allow for the possibility of outliers. This is achieved by introducing slack variables

Given the training data set

               
           
  0 i    for       T ii by     wx and   0 j   for       T jj yb     wx
. Thus, the optimization problem in (3.35) can be rewritten as:

                          2 1 T T 1 minmize 2 subject to , 0, 1, , n ii i i i i i i i ii C yb by in                               w wx wx (3.36)
The minimization in Equation (3.36) is a trade-off between the flatness of   f x (model complexity) and the degree to which deviations larger than  are tolerated. This trade-off is controlled by the user defined constant 0 C  . This method of tolerating errors is known as -insensitive loss function [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. The linear -insensitive loss function is described by Equation (3.37),, which is shown in the right part of Figure 3.6. Points that lie inside the - L is given as

                                        2 11 T 1 T 1 1 2 nn i i i i i i ii n i i i i i n i i i i i LC yb yb                                              w wx wx (3.38)
The Lagrangian L must be minimized with respect to ,

w , b   i  and   i   (the primal variables) and maximized with respect to   , i    , i     i  and   i   (the dual variables), where         , , , 0
i i i i       .
Then, the minimization of L with respect to the primal variables ( ,

w , b   , i    i  
) and the maximization with respect to dual variables involves finding the saddle point, at which the derivatives of L with respect to the primal variables have to vanish, i.e., 3.3 Surrogate Models 66

          1 0 n i i i i L               wx w (3.39)       1 0 n ii i L b          (3.40)       0 ii i L C         (3.41)       0 ii i L C           (3.
                                              ,1 11 1 1
,, 2 maximize subject to 0 and , 0,

n i i j j i j ij nn i i i i i ii n i i i i i Lk y C                                xx (3.43)
where

                      T ,, i j i j i j k        x x x x x x
is the kernel function. This maximize problem usually solved by quadratic programming algorithms, and thus   i  and

  i   are obtained.
From Equation (3.39) we can obtain

          1 n i i i i          wx (3.44)
Thus, the prediction of SVR for any new input point is computing by

            1 ˆ, n i i i i f k b        x x x (3.45)
where

                T ,, i i i k        x x x x x
x is the kernel function between the new input x and sampled point   i x in training data set. The constant term b , known as the bias, is obtained by exploiting the so-called Karush-Kuhn-Tucker (KKT) conditions [START_REF] Smola | A tutorial on support vector regression[END_REF][START_REF] Boyd | Convex optimization[END_REF].

The use of kernel functions (or kernels)   , k  make SVR models have the capability of capturing complicated landscape. The kernel is related to the transformation    x , which maps the input space to the so-called feature space (i.e.   :

 x ), by               ,, i j i j k   x x x x
. In other words, the kernels are the inner products (also known as dot products) of input points. Conventionally, kernels need to satisfies continuous, symmetric and positive definite conditions. Popular choices for   , k  are:

 Linear:           ,, i j i j k  x x x x
 p degree homogeneous polynomial:

            ,, p i j i j k  x x x x
 p degree inhomogeneous polynomial:

            ,, p i j i j kc  x x x x
 Gaussian:

            2 2
, exp

i j i j k     x x x x

Model Validation/Model Performance Assessment

The accuracy of a surrogate model is affected by the type of surrogate model and the quality and quantity of the dataset from which it is constructed [START_REF] Ryberg | A metamodel-based multidisciplinary design optimization process for automotive structures[END_REF]. Before we use the built surrogate models, the performance of the surrogate model need to be assessed and validated.

This section discusses several measures and methods to assess the accuracy of a surrogate model.

Fitting Error and Prediction Error

The most popular and simple way to evaluate the accuracy of a surrogate model is to examine its residual errors, i.e., the difference between the output from original simulation

              11 , , , , nn yy  xx .
The residual error thus can be expressed as

    f f  xx .
Residual error evaluated on the training points is referred to as the fitting error, while residual error at the set of points which are randomly generated and not be touched during the surrogate model constructing stages is called prediction error. In the following, several popular measures of model accuracy that are based on residual error and thus can be evaluated on whether fitting error or prediction error are presented. The first two measures are the root mean squared error (RMSE) and the maximum absolute error (MAE) [START_REF] Wang | Review of Metamodeling Techniques in Support of Engineering Design Optimization[END_REF], separately defined by

      2 1 1 m ii i RMSE y y m    (3.46)     max , 1, 2, , ii MAE y y i m    (3.47)
where m is the number of validation points,   i y is output from experiment, and   ˆi y is the output predicted by surrogate model. When RMSE and MAE are calculated on fitting error, the training points are used as validation points. While, a set of untried points are used when we assess the surrogate model by prediction error.

An additional error measure based on the residual is the coefficient of determination 2 R , which is defined by For surrogate models that interpolate the training points, there are no residual error on training points, i.e., the fitting error is zero. Consequently, measures based on fitting error are not relevant for interpolating surrogate models. In this case, the prediction error or the cross validation described in next subsection can be used to assess the performance of the model. Usually, prediction error, which assesses the ability of the surrogate model to predict responses in unknown design points, is a more useful measure than fitting error.

                 

Cross Validation

In many situations, computer experiments are computationally expensive. Thus, in such situations, evaluation of the prediction error is expensive or time-consuming. To alleviate this problem, a general strategy is to estimate the prediction error of a surrogate model using the cross validation (CV) procedure [START_REF] Fang | Design And Modeling for Computer Experiments[END_REF].

In cross validation (CV), the sample data is divided into training and test points. The training points are used to build the surrogate model, while the test points are used to test the performance of the model. The cross-validation technique operates through the following five steps [START_REF] Messac | Optimization in Practice with MATLAB®: For Engineering Students and Professionals[END_REF]:

1) Splits the sample points randomly into p (approximately) equal subsets;

2) Removes each of these subsets in turn (one at a time);

3) Trains a surrogate model according to the remaining 1 p  subsets; 4) Computes the error of the built surrogate model using the omitted subset; 5) Once each one of the p subsets has been used as the omitted subset, the p sets of errors are generally aggregated to yield a global error measure. Above described CV technique, which splits the sampling data into p subsets, is known as p-fold CV. A variation of p-fold CV is the leave-k-out approach, in which all possible subsets of size k are left out, and the surrogate model is constructed to the remaining set.

Each time, the error measure is evaluated at the omitted points. If 1 k  , the cross validation in this special case is called leave-one-out CV. Leave-one-out CV is probably the simplest and most widely used method for surrogate model validation when additional validation points are not available. The root mean squares prediction error for the one-leave-out CV is calculated by:

      2 1 1 n ii CV i i RMSE y y n     (3.49)
where   i y is the response at   i x from sample data and   ˆi i y  is the prediction at   i x from the surrogate model constructed by using the sample points except

      , ii y x
. The leave-one-out CV is a measure of how sensitive the surrogate model is to lost information at its data points.

However, an insensitive surrogate model is not necessarily accurate and an accurate model is not necessarily insensitive to lost information. Therefore, the leave-one-out CV is not

A Brief Survey of Surrogate-Assisted Evolution Strategies

The

Surrogate-Assisted Evolution Strategies

The first surrogate-assisted ES, very likely, were proposed by M. Papadrakakis in [START_REF] Papadrakakis | Structural optimization using evolution strategies and neural networks[END_REF],

where a hidden-layer ANN was used to predict the objective and constraint functions values of an expensive structural optimization problem using   In [START_REF] Bajer | Investigation of Gaussian Processes and Random Forests as Surrogate Models for Evolutionary Black-Box Optimization[END_REF], generation-based evolution control was adopted in surrogate-assisted CMA-ES,

where Gaussian process and random forests models were used as surrogates.

Pre The basic idea behind this approach is that only the most promising individuals with a good fitness prediction are evaluated with the true fitness function, which results in a reduction of the number of expensive true fitness calls. The Model-Assisted Evolution Strategy (MAES) using pre-selection strategy for   , -ES  was described in [START_REF] Ulmer | Evolution strategies assisted by Gaussian processes with improved preselection criterion[END_REF].

Applications of Kriging model based pre-selection strategy in ES were studied in [START_REF] Ulmer | Optimization by Gaussian Processes assisted Evolution Strategies[END_REF][START_REF] Ulmer | Evolution strategies assisted by Gaussian processes with improved preselection criterion[END_REF][START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF][START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF]. The authors suggested to use criteria that based on both the model prediction   f x and estimated standard deviation   ŝ x to identify the most promising individuals in preselection. This is based on the opinion that the key issue of using approximate models for evolutionary computation lies in the trade-off between the exploitation of the approximate model by sampling where it is optimized and the need to improve the approximate model by sampling where the model confidence is low [START_REF] Ulmer | Evolution strategies assisted by Gaussian processes with improved preselection criterion[END_REF][START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF].

In 2004, Runarsson [START_REF] Runarsson | Constrained Evolutionary Optimization by Approximate Ranking and Surrogate Models[END_REF] proposed the so-called approximate ranking procedure to evaluate the quality of the surrogate models. The approximate ranking procedure, which is described in Algorithm 4.

 

1:

1 i i     x where       1: 2: : ˆˆf f f        x x
  ˆ, 1, , kk f f k   x 8: determine the parent set   : 1 mi i     x where       1: 2: : ˆˆf f f        x x x 9: if 1 mm  
  1 , , , kk k tf   x
Recently, comparison-based surrogate models [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF], such as ordinal regression [START_REF] Runarsson | Ordinal Regression in Evolutionary Computation[END_REF] and ranking Support Vector Machine [START_REF] Loshchilov | Self-adaptive Surrogate-assisted Covariance Matrix Adaptation Evolution Strategy[END_REF], are used in surrogate-assisted ES. The rank-based Support Vector Machines (SVM) was used to learn the surrogate model for CMA-ES in [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF][START_REF] Runarsson | Ordinal Regression in Evolutionary Computation[END_REF][START_REF] Loshchilov | Self-adaptive Surrogate-assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]. The experimental validation demonstrated the invariance of monotonous transformation of the fitness by using rank-based SVM. However, the main weakness of comparison-based surrogates assisted ES is that it can not handle well multi-modal diversity, i.e., comparison-based surrogate models easily fail in accounting for multi-modal landscapes of fitness functions [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF][START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF]. evaluations ago [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF]. Thus, in this case, the built model would be very precise for the smaller cluster near the local optimum, but may not perform well in other place. This is undesirable in surrogate-assisted evolutionary search. When surrogates are built for each individuals (such as in [START_REF] Lu | An evolution strategy assisted by an ensemble of local Gaussian process models[END_REF]), usually the k-nearest neighbor points of each individual are used to train a surrogate model for corresponding individual.

Some Open Issues in Surrogate-Assisted ES

k-Nearest

Model Quality Assessment

Without no doubt, the most popular measure for model quality is the Mean Square Error (MSE) which assesses the mean squared difference between the individual's original fitness function   f x and the output of the surrogate model  

f x , that is       2 1 1 MSE . n ii i ff n    xx (4.1)
The MSE is evaluated over n different individuals that are taken into account for the estimation of the model quality, for instance, the n   offspring individuals in one generation.

The MSE is widely used to evaluate the approximation accuracy of surrogate models.

In the investigation of quality measures for surrogate model in evolutionary computation by Jin et al. [START_REF] Jin | Quality measures for approximate models in evolutionary computation[END_REF], it has been stated that the Mean Square Error (MSE) of the model only weakly where    0   is the number of correctly selected individuals, i.e., the number of individuals that have would also been selected if the original fitness function had been used in fitness evaluation. The expectation of in case of random selection, that is

2 0 , m mm m                            (4.
  max 1 1 . 2 m m                 (4.

Issues in Model Exploitation

The pre-selection and evolution control for model exploitation have already been mentioned in Section 4.  brings about a larger selection pressure and results in a stronger impact of the model f on the convergence of the optimization process [START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF]. In [START_REF] Ulmer | Optimization by Gaussian Processes assisted Evolution Strategies[END_REF][START_REF] Ulmer | Evolution strategies assisted by Gaussian processes with improved preselection criterion[END_REF][START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF][START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF], the pre-selection population size Pre  was set as a constant during the optimization. The concept of model impact control was proposed by Ulmer et al. [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF], where the value of Pre  , respectively. Commonly, in evolution control, the number  is set by the users. For example, in [START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF] generation-based evolution control was carried out every third generation. Empirical investigation on individual-based evolution control has performed in [START_REF] Jin | On evolutionary optimization with approximate fitness functions[END_REF], which showed that more than half of the individuals in population need to be controlled when the random selection strategy is used and about 40%

of the individuals should be controlled when using best strategy, in order to guarantee a correct convergence. To manage the control frequency of generation-based control, Jin et al.

        , for 1, , g g g kd k     x m 0 C
. Obviously, it can also be stated that the  offspring are from a multivariate normal (MVN) distribution with mean vector

  g  m  and covariance matrix       2 gg   C  , i.e.

 

, kd x   . Correspondingly, the Mahalanobis distance from a point x to current distribution mean

  g  m  is calculated by                       1 T2 T 1 ,. g g g g           x x x x m C x m     (4.7)
Generally speaking, data points with smaller Mahalanobis distance to current mean of mutation distribution are more relevant to current offspring population. Thus, we can select the k-Nearest Neighbor Points according to the Mahalanobis distance to current distribution mean   g m as the training set.

Confidence Interval for Training Set Selection

Considering the  offspring are from a MVN distribution with mean vector

  g  m  and covariance matrix       2 gg   C  , i.e.
                     
    max selection 1 1 . 2 i Qi            (4.10)
The expectation value of selection Q for a purely random selection process is given as the product of the expectation value of the number of correctly selected individuals and the expectation value of the rank of a correctly selected individual :

  2 rand selection 11 1 2 1 . 2 ii ii Q i i                                    
Q                  (4.12)
where Pre  is the adaptation rate. It was demonstrated in [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF] that controlling the impact of model impact ( Pre  ) enhances the performance of model assisted Evolution Strategies (MAES) with fixed model influence (fixed Pre  ).

Proposed Pre-Selection with Model Impact Control Proposed Selection-based Model Quality Measure using Recombination Weights

In evolutionary algorithms, the selection operator choosing the best individuals to enter next evolution loop and thus gives the evolutionary search a direction. Specially, in evolution strategies, the selection is based on fitness ranks and deterministic. The recombination after selection is also deterministic in ES. Thus, selection is particularly important in evolution strategies. Consequently, in surrogate-assisted evolutionary computation, the quality of model-based selection (selection by using surrogate model) is of importance, particularly in surrogate-assisted evolution strategies.

In the investigation of quality measures for surrogate model in evolutionary computation by Jin et al. [START_REF] Jin | Quality measures for approximate models in evolutionary computation[END_REF], it has been stated that the the Mean Square Error (MSE) of the model only weakly correlates with the ability to correctly select individuals in evolutionary computation, and selection based or ranking based model quality measurements were suggested to use. In [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF], a selection based model quality was used to control the size of pre-selection population for surrogate-assisted evolution strategies. This method for updating Pre  has pointed out a direction for pre-selection population size control, :

g i i i w      
mx , which is the new mean vector of mutation distribution. It apparent that individuals with high fitness ranking have larger weights and thus contribute more to the search step.

Considering the model based selection process selecting  out of the  individuals with the best predicted fitness, an individual is correctly selected, if it would also be selected by a selection process based on the original (true) fitness of the individuals. Usually, the number of correctly selected individuals can be used to measure the quality of the surrogate model.

We define a selection based model quality that takes in account not only the number of correctly selected individuals but also the real fitness ranking of individuals. Our idea is that, for each model-based selected individual (individual that is selected according to their prediction fitness), if it is also selected by selection process using real fitness, its associated

wi ii ii Q i w                                             (4.14)
In surrogate-assisted evolution strategies using pre-selection strategy, to measure the quality of the pre-selection process with above defined w Q , we need to know the true fitness of all Pre  individuals. However, only the true fitness of  most promising pre-selected individuals is known. We also taken the assumption in [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF] that the model quality w Q for a selection of  out of the  individuals, for which we know the predicted and the true fitness, is equivalent to that of the pre-selection process.

Controlling the Size of Pre-selection Population (Model Impact) 𝝀 𝐏𝐫𝐞

The evolutionary search can benefit from the surrogate model assistance if the model selection process performs better than a purely random selection and correspondingly model selection quality is better than the quality of a purely random selection.

After each generation in the surrogate-assisted ES with pre-selection, the actual measured selection quality

              (4.16)
The critical values of model quality measures are listed in Table 4.1. For model quality measure rank  , rank 0.5   indicates a moderate monotonic relationship between the predicted fitness and exact fitness. Thus, we choose rank 0.5

 

as the threshold value of rank  .

Maximum value max Q Threshold value TS Q w Q 0 1  selection Q 0   12         2 2 1 2        rank  -1 1 0.5

Algorithm of Pre-Selection Strategies

The pre-selection concept has been introduced in Section 4.

Individual-based Evolution Control

In surrogate-assisted evolution strategies with evolution control, the control frequency, which is described in Section 4.2.2.3, is of significance. On one hand, in order to prevent the evolutionary search from being misled by a false optimum introduced by the surrogate model [START_REF] Jin | A framework for evolutionary optimization with approximate fitness functions[END_REF], the surrogate model should be used together with the original fitness function. On other hand, the original fitness function should be used sparsely such that computational cost on (original) fitness evaluation is saved as much as possible. Evolution control could be regarded as a scheme to make balance between saving the computational cost of fitness evaluation and ensuring the evolutionary search converging to the global optimum or a nearoptimum of the original fitness function. This subsection and next subsection devote to putting forward strategies for individual-based and generation-based evolution control, respectively.

The concepts of evolution control have already introduced in Section 1.3.3. In surrogate-assisted evolutionary computation, the evolution control means that the original fitness function is used to evaluate some/all individuals in some/all generations. An individual that is evaluated using original fitness function is called a controlled individual.

Similarly, a generation in which all its individuals are evaluated using the original fitness owing to the adaptive control frequency, and thus may save more computational cost on fitness evaluation and consequently would be more efficient.

Fixed Individual-based Control (FIC) using Metric

In

    Mean Ĉ f  xx (4.17)     SD Ĉ s  xx (4.18)       SLB ˆĈ f A s    x x x (4.19)           min POI ˆîf 0 ˆ0 if 0 ff s s C s                x x x x x (4.20)                     min min min EI ˆˆ, if 0 ˆˆ0, if 0 f f f f f f s s ss C s                                xx x x x xx x x (4.21)
where

A is a constant and in this work we set 2 , , , , , , The original approximate ranking procedure proposed by Runarsson [START_REF] Runarsson | Constrained Evolutionary Optimization by Approximate Ranking and Surrogate Models[END_REF] has been presented in Algorithm 4.1. In the iteration of the original approximate ranking procedure, one individual is selected to be re-evaluated by original fitness function, and then the surrogate model is updated and the parent set is selected based on the updated surrogate.

A  ,    
given:             1 , , , , , , , ,
            min max 1 , , , , , , , , , , ggg kk
                    x x x x x //
This loop is stop until the parent set does not change or all the individuals in the generation have been re-evaluated. However, when the population size  is large, which is often required in solving difficult multimodal problems or high dimensional problems, the amount of information added in one iteration may result in insignificant changes even of a surrogate model with bad ranking predictions [START_REF] Kern | Local Meta-models for Optimization Using Evolution Strategies[END_REF]. To overcome this deficiency, we make some modifications for the original approximate ranking procedure.

Algorithm 4.5 Modified Approximate Ranking Procedure (ARP)

1: given:             1 , , , , , , , ,
ggg kk k t f C    z x m C x x 2: approximate:   T training_set_selection  ,   T ˆmodel_training f  and predict   ˆ, 1, , kk f f k   x
3: rank and determine the parent set

 

1: proportional to  , which is similar to that in [START_REF] Kern | Local Meta-models for Optimization Using Evolution Strategies[END_REF]. Furthermore, since Kriging model is used, the metric (including Mean, SD, SLB, POI and EI) can be adopted in selecting individuals for re-evaluation and updating the model. With the aid of metric, the individuals with which the model can be significantly improved can be selected for re-evaluation and added to the training data set. Therefore, with above modifications, the efficiency of approximate ranking procedure could be enhanced. The modified approximate ranking procedure is expressed by Algorithm 4.5.

1 i i     x where       1: 2: : ˆˆf f f        x x

Generation-based Evolution Control

For generation-based evolution control, both fixed and adaptive control frequency can be used. In fixed generation-based evolution control, the control frequency is fixed, i.e., generation control is carried out once in a fixed number of generations. While, for adaptive generation-based evolution control, the frequency of generation control is adaptive, usually depend on the fidelity (quality) of the surrogate model. An appropriate generation control frequency can ensure the evolutionary computation converge to correct optimum and, at the same time, reduce the computational cost as much as possible.

In [START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF], generation control was performed in every third generation in the comparative studies of surrogate-assisted ES with pre-selection strategy and generation-based control. Jin et al. [START_REF] Jin | Managing approximate models in evolutionary aerodynamic design optimization[END_REF][START_REF] Jin | A framework for evolutionary optimization with approximate fitness functions[END_REF] proposed a strategy for adapting the generation control frequency based on the fidelity of the surrogate model. This adapting strategy for control frequency is presented below.

It is intuitive that the higher the fidelity of the surrogate model is, the more often the fitness evaluation can be made using the surrogate model [START_REF] Jin | Managing approximate models in evolutionary aerodynamic design optimization[END_REF]

        2 1 1 îi i E k f f      xx (4.23)
where  is the population size of ES, , 1, 2, , i i

  x
are  individuals that are evaluated using original fitness function in current control cycle. Apparently,  

Ek is the mean square error of  evaluated points. To make sure that the information on the model fidelity is always available, at least one generation should be controlled within one evolution control cycle. Besides, since the model fidelity is estimated locally based on the error information from the last cycle, p should not be too large.

Proposed Adaptive Generation-based Control (FGC) based on Model Quality

In previous sections, the model quality has been used to control the model impact of pre-selection strategy. In this section, it is proposed to use the model quality to determine

whether

Variable Transformation for Model Learning and Prediction

Recently, the covariance matrix information has been used in model training so that the correlations between variables are taken into account during model learning. In lmm-CMS-ES (local-meta-model CMA-ES) [START_REF] Kern | Local Meta-models for Optimization Using Evolution Strategies[END_REF] and its relevant studies [START_REF] Bouzarkouna | Investigating the Local-Meta-Model CMA-ES for Large Population Sizes[END_REF][START_REF] Kern | Fast quadratic local meta-models for evolutionary optimization of anguilliform swimmers[END_REF], the covariance matrix C that is adapted in CMA-ES was used as a metric in the calculation of distance between

two variable vectors       T 1 , i j i j i j d     x x x x C x x ,
        , , 1, , g g g kd k     x m 0 C
. In other words, the  offspring are from MVN distribution with mean vector

  g  m  and covariance matrix       2 gg   C  , i.e.,             2 ,, g g g k d d  x m C  
. Thus, for the original variable vector

d i  x , the following transformation T is used               12 2 12 T, g g g i i i i            x x x C x m   (4.24)
where

d i   x
is the associated vector of

d i 
x obtained by the defined transformation. This transformation uses the inverse of the covariance matrix of MVN distribution and has the effects of ( 1) standardizing all variables in the vector to the same variance and (2) eliminating correlations between variables [START_REF] Rencher | Methods of multivariate analysis[END_REF]. It is worth mentioning that standardization of data set is a common requirement for many machine learning implementations. With transformation in Equation (4.24), the Mahalanobis distance in Equation (4.7) can be rewritten as

              T TT 1 , T T . i i i i         x x x x x x x     (4.25)
This shows that the Mahalanobis distance between original variable vectors can be viewed as the Euclidean distance on transformed vectors.

Initial Sampling and Informed Start Point

For surrogate-assisted ES, a set of evaluated points is indispensable before the starting of surrogate assistance. Usually, the surrogate model is trained at the beginning with a randomly created initial population (or initial sampling), which can be found in [START_REF] Ulmer | Optimization by Gaussian Processes assisted Evolution Strategies[END_REF][START_REF] Ulmer | Evolution strategies assisted by Gaussian processes with improved preselection criterion[END_REF][START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF][START_REF] Ulmer | Model Assisted Evolution Strategies[END_REF].

In this work, a different way is used to create initial population for constructing the initial surrogate mode. As stated in Chapter 3, space-filling design is appropriate to represent all portions of the design space. Thus, in this work, the initial sampling is created by the spacefilling Latin hypercube designs described in Section 3.2.2 rather than by randomly sampling. It is apparent that smaller number of iterations are needed to reach the optimum when the start point   0 m for the evolution loop is near the location of the optimum, and vice versa.

However, for black-box problem, usually no prior knowledge about the global optimum is available before we solve the optimization problems. Thus, the start point

  0 d  m is
generally initialized uniformly and randomly within the search space. Fortunately, information can be extracted from the samples of initial sampling (DOE). Therefore, we proposed to choose the best point in initial sampling (the point with minimum fitness function value) as the start point for subsequent evolution loop, after the evaluation of the DOE n points in DOE. This strategy for start point selection is referred to as informed start point in this work. In this work, the initial sampling and informed start point are performed before the evolution loop for all surrogate-assisted evolution algorithms.

KA-CMA-ES using Pre-Selection

This subsection dedicates to the KA-CMA-ES using pre-selection, where pre-selection

2: initialize

        0 0 0 0 cT 0, , , , = , , 0, 0 gt            p 0 p 0 C I 3: initial sampling:   DOE T T DO D 1 E E O lhsdes , ign , , n nd    X x x
4: evaluate initial samples:

    ,, k k k k f f f  xx where DOE DOE , 1

8:

model training:

  T ˆmodel_training f  9:
if PS is used then 10:

          Pre 1 , , pre-selection_procedure , , , , ggg k k k k f     x z m C //pre-selection without model impact control 11: elseif CPS is used then 12:             Pre 1 , , pre-selection_procedure , , , ,
g g g g k k k k f     x z m C //
      1 1 model_qualtiy_estimate , , ,
g k k k k k k Q f f         xx 21:             TS TS Pre Pre max TS 1 Pre TS Pre Pre TS min , if , otherwise g gg g g g QQ QQ QQ QQ QQ                   //update Pre

KA-CMA-ES using Individual-based Evolution Control

2: initialize

        0 0 0 0 cT 0, , , , = , , 0, 0 gt            p 0 p 0 C I 3: initial sampling:   DOE T T DO D 1 E E O lhsdes , ign , , n nd    X x x
4: evaluate initial samples: ,,, , , fixed_individual_control , , , , , nn and metric  

    ,, k k k k f f f 
              1 1 , , , ,
ggg kk k k k k k tf t f C            z x m C xz xx 15: elseif MIC is used then 16:               1 1 min max , , , , , , , , , mixed_individual_control , , , , ,
ggg kk k k k k k tf t f C            z x m C xz xx 17: end if 18: update      
C x ). 2: initialize         0 0 0 0 cT 0, , , , = , , 0, 0 gt            p 0 p 0 C I 3: initial sampling:   DOE T T DO D 1 E E O lhsdes , ign , , n nd    X x x
4: evaluate initial samples:

    ,, k k k k f f f 
                11 , , , , approximate_ranking_procedure , , , , , , , ,
ggg k k k k k kk t f t f C     x z z x m C x x 12: update      
: initialize           0 0 0 0 0 c T control 0, , , , = , , 1, 0, 0 I g t             p 0 p 0 C I 3: initial sampling:   DOE T T DO D 1 E E O lhsdes , ign , , n nd    X x x
4: evaluate initial samples:

    ,, k k k k f f f 
      1 1 model_qualtiy_estimate , , , g k k k k k k Q f f         xx 17: if   TS g QQ  then 18:   1 control 0 g I  
//mark next generation as model-evaluated 19: else 20:

  1 control 1 g I  
//mark next generation as controlled 21:

end if 22: esleif   control g I is equal 0 then 23:   ˆ, 1, , kk f f k   x

24:

 

Experimental Studies

The validation and the performance evaluation of optimization algorithms is commonly carried out by using a chosen set of benchmarks or test functions. In this section, the proposed Kriging-assisted CMA-ES algorithms in Section 4.4 are validated using a set of test functions and their performance are evaluated and analyzed.

Experimental Setup

Test Functions

Test functions are important in performance validation and comparison of optimization algorithms. In order to comprehensively evaluate an algorithm, the used set of test functions should includes enough functions with different characteristics, such as continuous, discontinuous, unimodal, multi-modal, separable, non-separable. In this work, we focus on continuous optimization problems. In the experimental studies of proposed Kriging-Assisted CMA-ES (KA-CMA-ES) algorithms, a set of 12 continuous benchmark functions is carefully selected from [START_REF] Jamil | A literature survey of benchmark functions for global optimisation problems[END_REF][START_REF] Suganthan | Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization[END_REF][START_REF] Chen | Problem definitions and evaluation criteria for CEC 2015 special session on bound constrained single-objective computationally expensive numerical optimization[END_REF].

All the test functions are minimization problems defined as follows : . For handling the box constraints, the resampling method, i.e. re-sampling any infeasible solution x until it becomes feasible, is adopted.

      T 12 LB UB minimize , , ,
The name, expression and the search space of 12 test functions are listed in

D i i fx    x   5, 5 D  Bent Cigar   2 4 2 21 2
10

D i i f x x    x   100,100 D  Sum Squares     2 3 1 D i i f i x    x   10,10 D  Schwefel 1.2   2 4 11 Di j ij fx       x   100,100 D  Powell Sum   1 5 1
D i i i fx     x   1, 1 D  Schwefel Absolute   6 1 1 D D ii i i f x x      x   100,100 D  Rosenbrock       1 2 2 2 71 1 100 1 D i i i i f x x x            x   5, 5 D  Ackley       2 8 11 11 20exp 0.2 exp cos 2 20 exp 1 DD ii ii f x x DD                 x   32, 32 D  Levy             1 2 22 91 1 2 2 sin 1 1 10sin 1 1 1 sin 2 D ii i DD f z z z zz                  x   10,10 D  Weierstrass         max max 10 1 0 0 cos 2 0.5 cos 2 0.5 D k k k k k k i i k k f a b x D a b                       x where 0.5, 3, max 20 a b k      0.5, 0.5 D  Bohachevsky       1 22 11 1 1 1 2 0.3cos 3 0.4cos 4 0.7 D i i i i i f x x x x             x   15,15 D  Rastrigin     2 12 1 10 10cos 2 D ii i f D x x         x   5, 5 D 

Experimental Setting

For each test problem, 25 independent runs are performed using the algorithm which are chosen to be investigated. The initial mean is set as: for CMA-ES, the initial mean   0 m is uniformly sampled within the search domain of the problem (given in Table 4.2) ; for KA- In experimental study, the performance of an algorithm is evaluated according to the success rate (SR), the success performance (SP) and the speedup performance (SPU). The speedup performance (SPU) denotes the degree of improvement in performance (according to number of exact function evaluations) of the KA-CMA-ES algorithm.

Obviously, the SPU of an algorithms is larger than one (SPU>1) means the algorithm makes improvement in success performance. On the contrary, SPU<1 indicates the algorithm performs worse than standard CMA-ES. From the radar chart of average success rate (SR) in left of Figure 4.3 and Table 4.4, for each category of problems, the KA-CMA-ES using PS have average success rate exceeds 0.9. Specifically, the PS-kNN has the highest average success rate on problems with D=2 and D=20, and unimodal problems. PS-Interval has the highest average success rate on problems with D=5, D=10, multimodal problems, and on overall problems. The PS-Recently method has the lowest average success rate among three methods. PS-kNN and PS-Interval has higher average success rate than the standard CMA-ES on all categories of problems. In terms of speedup performance (SPU) or success performance, PS-Interval has the best success performance on all other category problems except 2 dimensional problems. Thus, it can also be found that the 'Interval' method is superior to 'Recently' and 'kNN' methods. The experiments results of KA-CMA-ES using pre-selection with model impact control (CPS), i.e., the success rate (SR), success performance (SP) and speedup performance (SPU), are listed in Table 4.5. In the table, CPS-

Experiments on KA-CMA-ES using

Comparison of Pre-Selection with and without Model Impact Control

The results of KA-CMA-ES using pre-selection without model impact control using 'Interval' training set selection method (PS-Interval) and KA-CMA-ES using pre-selection with model impact control are summarized in Table 4.7. It is apparent that pre-selection with model impact control (CPS) outperforms pre-selection without model impact control (PS)

for most of the test problems according to the speedup performance. Thus, model impact control is suggested in KA-CMA-ES using pre-selection.

The average success rate and speedup performance of PS-Interval, CPS- The results of KA-CMA-ES using Fixed Individual-based Control (FIC) with different metrics (Mean, SD, SLB, POI and EI) are given in Table 4.9. In all the tables of experiments results in this chapter, we mark the SPU that do not exceed one ( SPU 1  ) in red bold format, which means the corresponding algorithm does not improve the success performance on the corresponding test problem, and mark the maximum SPU of each row in bold number, which represents the best success performance among the investigated algorithms on the corresponding problem. Additionally, the average success rate and speedup performance of FIC using five different metrics on each category problems and overall problems are presented in Table 4.10 and plotted in Figure 4.6. In the tables, FIC-Mean stands for fixed individual-based control using Mean as metric, FIC-SD denotes FIC using SD as metric, FIC-SLB means using SLB as metric in fixed individual-based control, FIC-POI signifies the FIC using POI metric, and FIC-EI represents fixed individual-based control using EI as metric.

It is no doubt that the FIC-EI, on the whole, has the highest success rate among the algorithms of FIC using five metrics. With regards to the success performance and speedup performance, FIC-EI also has the best performance. Among the five metric, these metric (SLB, POI and EI) which balances the exploitation of the surrogate and exploration of search space performs better than metrics that only consider exploitation or exploration (Mean, and SD). Taking both SR and SP into consideration, FIC-EI has advantages over other metrics and thus is suggested to be used in other KA-CMA-ES algorithms using metric, such as KA-CMA-ES using approximate ranking procedure which will be investigated in section 4.5.4.

Mixed Individual-based Control (MIC)

In the proposed mixed individual-based control (MIC), two metrics, Mean and SLB, can be used. The results of KA-CMA-ES using MIC are given in Table 4.11. The average success rate and speedup performance on each category of problems are listed in Table 4.12 and shown in In Table 4.11, in terms of speedup performance, MIC-Mean (mixed individual-based control using metric Mean) performs worsen than standard CMA-ES on 10 dimensional 6 f and 10 f and 2 dimensional 12 f . While, MIC-SLB (mixed individual-based control using metric SLB) has better success performance than standard CMA-ES on all the test problems.

Form this, it can be found that MIC-SLB is more stable and reliable than MIC-Mean.

From Figure 4.7, it is apparent that MIC-SLB has higher success rate (SR) than that of MIC-Mean, especially on multimodal problems. Furthermore, generally, MIC-SLB has better average speedup performance (SPU) than MIC-Mean on all categories of problem except problems with D=5. Taking both the success rate and speedup performance into consideration, it can be concluded that MIC-SLB outperforms MIC-Mean.

Comparison of Fixed and Mixed Individual-based Control

We compare the performance of KA-CMA-ES using Fixed Individual-based Control (FIC) and KA-CMA-ES using Mixed Individual-based Control (MIC). The success rate (SR) and speedup performance (SPU) of KA-CMA-ES using FIC and MIC are listed and compared in Table 4.13. The average SR and SPU of each category problems of KA-CMA-ES using FIC and MIC are given in Table 4.14 and shown in Figure 4.8.

From Table 4.13, it can be seen that, for most of the problems, KA-CMA-ES using MIC (the proposed mixed individual-based control) has higher speedup performance (better success performance) than using FIC (fixed individual-based control). It worth to note that for 10 dimensional 6 f , only FIC-EI and MIC-SLB perform better than the standard CMA-ES, and MIC-SLB also outperforms FIC-EI on this problem. Therefore, the proposed mixed individual-based control (MIC) has better success performance than fixed individual-based control (FIC).

Additionally, from the radar chart of average SR in the left of

Experiments on KA-CMA-ES using Approximate Ranking

Procedure

The experimental results of KA-CMA-ES using FIC with different metrics have proved that the metric EI is outstanding among the five metrics. Thus, in KA-CMA-ES using The results of KA-CMA-ES using ARP are listed in Table 4.15. The average success rate and speedup performance are given in Table 4.16 and plotted in Figure 4.9. In the tables, ARP-Mean indicates in Mean is used as metric in the ARP, and ARP-EI represents ARP using EI as metric.

From Table 4.15, we can find that ARP-EI (approximate ranking procedure using EI metric) outperforms that using ARP-Mean. On problem of 6 f with D=10, both ARP-Mean and ARP-EI perform worsen than CMA-ES.

According to the the radar chart of average success rate (SR), the average SR of ARP-EI are higher that of ARP-Mean on all categories of problems. In terms of the speedup performance, ARP-EI also performs better than ARP-Mean. Therefore, it can be concluded that ARP-EI outperforms ARP-Mean. The success rate (SR) and speedup performance of above selected KA-CMA-ES algorithms are listed in Table 4. [START_REF] Jin | Surrogate-assisted evolutionary computation: Recent advances and future challenges[END_REF]. Additionally, the average success rate (SR) and speedup performance (SPU) of each category of problems are given in Table 4.20 and illustrated by among the investigated algorithms of CPS, while rank CPS- is more stable and has higher success rate. With comparison of KA-CMA-ES using PS and CPS, it is apparent that KA-CMA-ES using CPS significantly perform better than that using PS.

In the experiments of fixed individual-based control (FIC), five different metrics have been investigated. The results demonstrated that the EI metric has the best performance among these five, considering both the success rate and success performance. For the proposed mixed individual-based control (MIC), the MIC using SLB metric outperforms that using Mean as metric. By comparing KA-CMA-ES using FIC and MIC, KA-CMA-ES using MIC significantly outperform FIC according to the speedup performance and success material modeling, currently, continuum damage mechanics (CDM) provides a comprehensive framework for the modeling of material's behavior with consideration of damage from a point view of continuum mechanics [START_REF] Murakami | Continuum Damage Mechanics[END_REF]. Since an aluminum material is examined in this paper, the ductile damage behavior of the material is modeled by the application of continuum damage mechanic theories.

Damage is defined as the presence and evolution of cracks and cavities at microscopic, mesoscopic or macroscopic level of materials which result deterioration in mechanical properties and may, eventually, lead to failure [START_REF] Souza Neto | Computational methods for plasticity: theory and applications[END_REF]. According to the difference in materials and loading conditions, damage may be generally divided into: ductile damage, brittle damage, creep damage, and fatigue damage. Ductile damage, occurring simultaneously with large plastic deformation, is the most typical damage in engineering fracture problems. So, in this work, ductile damage is considered and coupled into an elastic-plastic model.

In the theory of damage mechanics, the deterioration of ductile damaged material is generally assumed to be a process of voids nucleation, growth and coalescence. In the phase of voids nucleation, the micro-defects nucleate and germinate, i.e., the deterioration appears.

However, it has no significant effects on load-carrying capability of the material, and the damage effects of the material can be neglected. In the following phase of voids growth, the voids and cracks grow considerably as the plastic strain accumulating. The material properties are strongly affected in this phase and damage effects in mechanical properties cannot be ignored. During last phase of voids coalescence, the coalescence of voids and cracks arise and, eventually, the macroscopic crack if formed. The ductile fracture is induced at the end of this phase. The load-carrying capability of the material is reduced as the accumulation of ductile damage in the last two phases until the material losing its loading capability [93].

The continuum mechanics approach to ductile damage problems can be divided into the following two approaches: micromechanical approach and phenomenological approach.

In the first approach, the mechanical effect of damage is represented by the void volume fraction f . Then the plastic constitutive equation of voided material and the evolution law Non-gradient methods do not have limitations on objective function and their convergence speeds are enhanced significantly with the developments of non-gradient optimization techniques.

In the optimization problem of parameter identification using inverse method, the objective function is, normally, a highly nonlinear function and may possess many local minima [100]. Furthermore, the analytical form of the objective function is not known and the evaluation of the objective function is performed by numerical simulation (generally finite element method simulation), in other words, derivatives of the objective function are not available. Taking above aspects into account, a non-gradient global optimization algorithm is preferred to solve the optimization task in parameter identification.

As above mentioned that the objective function in parameter identification problem is not explicitly known and need to be evaluated by repeatedly running numerical simulations, this brings bout heavy computational burden and makes parameter identification process time-consuming. Therefore, in this chapter, the previously proposed KA-CMA-ES using ARP-EI, which has the best success performance among the investigated KA-CMA-ES algorithms, is applied in inverse method of parameter identification.

Elastic-Plastic Damage Model

In general, a normal elastic-plastic constitutive model contains the following basic components [START_REF] Souza Neto | Computational methods for plasticity: theory and applications[END_REF]:

(1) The elastic law, which gives the relation between stress and elastic strain;

(2) The yield criterion, which define the limit of plastic behavior;

(3) The flow rule, which describes the evolution of the plastic strain; and

(4) The hardening law, which characterizes the development of the criterion for subsequent yielding.

Accordingly, the elastic-plastic damage model, which considers damage in elasticplastic material, contains an evolution law of damage variables in addition to these above four components.

In this work, the commonly used von Mises yield criterion, the Swift strain hardening law and Lemaitre's ductile damage are included in the elastic-plastic damage model. The basic components of the elastic-plastic damage model are detailed in this section. In first place, the basic concept of damage variable is introduced.

Damage Variable

The notion of continuum damage mechanics was firstly proposed by Kachanov [101] when he studied on the brittle creep rupture of metal. The first damage variable with physical significance was given later by Rabotnov [START_REF] Souza Neto | Computational methods for plasticity: theory and applications[END_REF] who proposed the reduction of the crosssectional area due to micro-cracking as a suitable measure of the state of internal damage.

According to the concepts proposed by these authors, the damage variable, D , can be modeled by effective area reduction: From Equation (5.1), the damage variable D n is associated with n when cracks and cavities are oriented, which leads to tensorial nature of damage variable. However, if cracks and cavities are equally distributed in all directions, D n does not depend on n and it becomes a scalar D , which is called scalar damage variable. In this paper, the damage state is taken to be isotropic and the scalar damage variable D is used in damage model.

In the case that the state of damage is isotropic, the definition of damage by effective area reduction has allowed to define the so-called effective stress:

, 1 D     (5.2)
where is the effective stress, D is the scalar damage variable, and  is the second-order Cauchy stress tensor. The effective stress of Equation (5.2) simplifies the damage theory, and is employed in a number of damage problems, including ductile damage especially.

Constitutive Model

Lemaitre's ductile damage model is based on CDM theory and in the framework of thermodynamics of irreversible processes. The constitutive equations are derived from thermodynamic (or state) potentials and the evolution equations are given by dissipation potentials, which are detailed in Lemaitre's original model and CDM theory [START_REF] Murakami | Continuum Damage Mechanics[END_REF]. The elastic-plastic damage model used in this paper (with von Mises yield criterion, Swift hardening law and Lemaitre's ductile damage) are presented in the following.

To begin with, the strain tensor additive split is given: where  is the plastic multiplier, s is the deviatoric stress tensor. Accordingly, the evolution of accumulated (equivalent) plastic strain p  is:

.

1 p D    
(5.9)

The evolution equation of damage variable is denotes the trace, q is the von Mises equivalent stress and p is the hydrostatic stress.

The loading /unloading conditions 0, 0, 0       must be satisfied.

The equations of the elastic-plastic damage model in this chapter are summarized in Table 5.1. Loading/unloading conditions 0, 0, 0

     

Numerical Implementation of the Constitutive Model

The fully implicit elastic predictor/return-mapping scheme is commonly used in the numerical implementation of constitutive model [START_REF] Souza Neto | Computational methods for plasticity: theory and applications[END_REF]. Given the increment of strain The elastic predictor/return-mapping algorithm starts with the so-called elastic trial step.

In the elastic trial step, we assume that the step   1 , nn tt  is elastic; hence, neither damage nor strain hardening evolution takes place at this stage. Therefore, the elastic trial solution is

given by: The corresponding trial stress can be evaluated by the deviatoric/hydrostatic split of the stress tensor With Equation ( 14), the elastic trial von Mises equivalent stress can be computed:

trial trial trial 1 1 1 , n n n p    I  s (5.13) with             trial trial trial trial trial trial 1 1 d 1 1 1 1 trial trial trial trial trial trial 1 1 v 1 1 1 1 1 2 1 1 , 1 1 1 , e n n n n n n n e n n n n n n n D G D D p D K D p D p                           s s s ( 5
        trial trial trial trial trial trial 1 2 1 2 1 1 1 1 3 3 1 1 , n n n n n n q J J D D q            ss (5.16)
where trial

trial trial 3 1 1 1 2 : n n n q   
 ss is the effective elastic trial von Mises equivalent stress, which is evaluated from effective trial deviatoric stress. The yield function in the present case is then evaluated as In return-mapping step, the following nonlinear return-mapping equations need to solve:

    trial trial trial trial 1 11 trial 1 . 1 pp n y n n y n n q q D              (5.
          trial 1 11 11 1 1 1 11 1 1 1 1 3 =, 12 , 1 ˆ, 1 0. 1 ee n nn nn pp nn n s pp n n n n D n p n yn n D D Y D D H Dr q D                                                 s s (5.18)
The above system of nonlinear equations is unattractive for numerical implementation due to its high computational burden. Therefore, some relatively straightforward operations are performed to reduce the system of equations. According to the elastic law, together with Equation (5.19), the hydrostatic and deviatoric stress can be updated as

    1 1 1 trial 1 1 1 1 1 1, 3 1 2 , 2 n n n n n n n n p D p DG               s ss s (5.20)
where we have defined

  trial trial trial 1 1 1 1 1 1 trial trial 1 1 33 1 2 1 . 2 n n n n n n n n G D G D q                      s s s s s (5.22)
Thus, with the definition of von Mises equivalent stress, we obtain

  trial 1 1 1 1 3 . n n n q D q G         (5.23)
Consequently, the plastic consistency condition becomes

D q G p YD G D K                 (5.26)
Now, the return-mapping equation system (5.18) has been reduced to a system of two scalar equations (5.24) 5.2.

G fq DD YD f D D Dr                                     (5.
G fq DD YD f D D Dr                                       where     0 n pp y A            2 2 trial 11 1 11 2 1 13 , 6 (1 ) 2 nn n nn n D q G p YD G D K                 (iv) Update state         trial 1 1 1 trial 1 trial 1 1 1
The implementation of the fully implicit elastic predictor/return-mapping scheme is performed in ABAQUS TM subroutine VUMAT (a user subroutine to define material behavior) by FORTRAN codes. The purpose of parameter identification procedure is to obtain as good agreement as possible between simulated and experimental responses [105, 106]. In other words, the parameter identification process is to find a set of parameters which makes the minimum difference between simulation response and experiment data. Thus, inverse method is a useful tool for determining material parameters [107]. The inverse method of parameter's identification considers the parameter identification as an optimization problem, which is to find a set of parameters that minimize the difference between the experimental data and the 5.3 Inverse Method for Parameter Identification 162 numerical simulations. As an optimization problem, the inverse method of parameter identification basically consists two main parts: the first part is the formulation of the objective function,   f p , which measures the difference between the experimental data and numerical results; another part is the selection of an optimization strategy, which is able to find the minimum of the objective function.

Inverse Method for Parameter Identification

Framework of Inverse Method for Parameter Identification

The general framework of the inverse method for parameter identification is shown in Because the objective function is usually implicit and nonlinear, gradient-free global optimization algorithms are commonly used. However, the main drawback of the inverse method is its expensive cost in computation, which is due to the expensive evaluation of the objective function and large number of iterations before convergence. This has promoted the use of the KA-CMA-ES in inverse method.

Objective Function

In this work, our goal is to find a set of material parameters that yields the simulation force-displacement response with minimum difference between the experimental forcedisplacement curve, as shown in Figure 5.3. The function which provides a scalar measurement of the error between the experimental data and numerical simulation results is chosen as the objective function. We use a formulation based on a least square equation

[108], which is expressed as: In strain hardening parameter identification, the objective function   f p is computed only on the hardening part, i.e., the integral domain is from D0  to the necking point. In this situation, the softening part of force-displacement curve is not considered in objective function since damage is not included. When dealing with damage parameters, the softening part of the force-displacement curve needs to be taken into account. In this situation, we need to give precaution on the softening part, in order to build a right objective function.

                  2 sim exp 2 exp 2 sim exp 1 1 2 exp 1 1 F F D FD F F D D F D D n i i i i i n i i i i d f d                    p p p ( 5
In the numerical investigation of ductile damage parameters identification by Roux et

Parameter Identification using Inverse Method

The inverse method presented in the previous section is employed to identify the material parameters that are used in the ductile damage model presented in Section 5.2. The KA-CMA-ES using ARP-EI, which has been proven to be the outstanding algorithm among KA-CMA-ES that have been investigated in previous chapter, is adopted as the optimization algorithm in inverse method of parameter identification. At the same time, the standard CMA-ES is also used, in order to validate the results of inverse method using KA-CMA-ES and evaluate its performance. The objective function, which measures the difference between simulation and experimental response and is minimized in parameter identification, is defined in Equation (5.29). The maximum computational budget is set as 600, i.e., the maximum number of exact objective function evaluations is 600, which is a stopping criterion of the optimization process. Another used termination criterion for optimization is Poisson's ratio are: E=72.4 GPa and v=0.33. The geometry dimensions of the specimens are illustrated in Figure 5.5. The force-displacement curves can be directly obtained from the test system. The first part of the curve is the elastic strain, followed by plastic strain, and then necking and force decrease till the final fracture. The objective of parameter identification is to find an appropriate set of material parameters that yields the best forcedisplacement response, which has minimum difference with the experimental forcedisplacement curve.

In parameters identification process, the tensile test process is firstly modeled by Finite Element Method (FEM) simulation using the software ABAQUS. The test specimen is modeled by 3D deformable solid with the element type C3D8R. In the central part, the element size is 1.5 mm, there two layers in the thickness of the specimen. The FEM model of the tensile specimen is shown in Figure 5.6. The loads are applied at the two sides of the

Strain Hardening Parameters Identification

The Swift hardening law, which is widely used to model the strain hardening of metallic material, is used to describe the strain hardening behavior of the material. The yield stress defined by Swift hardening law is expressed in Equation (5.7). Previously presented inverse method is used to identify the parameters, A , 0  and n , in Swift law.

The inverse parameter identification of strain hardening parameters are independently carried out with the standard CMA-ES and KA-CMA-ES using ARP-EI. In the identification process, FEM simulations of tensile test are repeatedly performed to evaluate the objective function. The lower and upper bounds for the parameters in the identification process and the identified strain hardening parameters are listed in Table 5. with the identified parameters accurately describe the strain hardening behavior of the material. Furthermore, from Table 5.3, it can be seen that the parameter identification results from the inverse method using CMA-ES and those from the inverse method using KA-CMA-ES (ARP-EI) are almost the same. This can also be seen in the iteration history of A , 0  and n in Figure 5.7, where these three parameters of the identification process using CMA-ES and KA-CMA-ES converge at almost the same values. In other words, the inverse method using KA-CMA-ES gives a reliable solution for parameter identification problems.

The numbers of objective function evaluations (FES) of the identification process using CMA-ES and KA-CMA-ES are presented in Table 5.3. The identification process using the standard CMA-ES requires 533 objective function evaluations (533 FEM simulations of tensile test with different parameters); whereas, only 211 function evaluations are required in identification using KA-CMA-ES. Obviously, the computational cost of inverse parameter identification using KA-CMA-ES is significantly lower than that of identification using the standard CMA-ES. This can be found in the convergence graph in Figure 5.7, too.

The identification using KA-CMA-ES converges faster than that using the standard CMA-ES. Therefore, efficiency of the parameter identification process has been enhanced by using the proposed KA-CMA-ES.

Damage Parameters Identification

After the strain hardening parameters are identified in previous subsection, we turn to the identification of ductile damage parameters, p D  , r and C D . In the ductile damage model, von Mises yield criterion and Swift hardening law are used. The above identified parameters of Swift hardening law are used here. The lower and upper bounds for the parameters in the identification process and the identified strain damage parameters are listed in Table 5.4.

The convergence graphs and iterative history of the damage parameters are illustrated in From the results of damage parameter identification in Table 5.4, the identified damage parameters from inverse method using CMA-ES and KA-CMA-ES are very closed. This additionally indicates the KA-CMA-ES gives reliable results. The iteration histories of damage parameters in Figure 5.9 also demonstrate this.

In Figure 5.10, the force-displacement curve from the numerical simulation using the presented elastic-plastic damage model with the identified parameters is consistent with that from the tensile test. Additionally, the simulation results of the elastic-plastic model and elastic-plastic damage model, and the experimental data are also plotted in Figure 5.11.

Obviously, since damage effect is not considered in the elastic-plastic model, the difference between the experimental response and the elastic-plastic model results is significant after the occurrence of necking. In conclusion, the presented elastic-plastic damage model is capable of, and appropriate for, modeling material behavior with consideration of strain hardening and damage effect, and the proposed KA-CMA-ES using ARP-EI is reliable and efficient in inverse method of parameter identification.

Summary

             ,         1 ln 2 1 2 ln for 1, , ln 2 1 2 ln i j i wi j              ,   2 2 11 w i i ii ww      , 2 5 w w c d       , 1 1 2max 0, 1 1 w dc d             , 4 42 w c w d c dd      ,   1 2 2 1.3 w c d    ,   1 2 21 min 1 ,2 . 2 ww u w cc d          2: initialiser           0 0 0 0 0 c , 0, , , , 0 d g

Troisième Chapitre: Modélisation de Substitution

Ce chapitre examine la mé thode de modé lisation substitutive, qui est une approche pour dé crire le comportement (entré e-sortie) du modè le de simulation. Les trois é tapes de la modé lisation de substitution sont dé crites dans ce chapitre : (i) la conception des expé riences (DOE) ; (ii) la formation du modè le de substitution et (iii) la validation du modè le. Avec un modè le subrogé formé , la sortie des points non testé s peut ê tre pré dite à moindre coût par le modè le. Par consé quent, les modè les de substitution peuvent ê tre utilisé s dans l'optimisation é volutive en remplaç ant l'é valuation coûteuse de la fonction de fitness par la pré diction du modè le de substitution moins coûteux que l'é valuation exacte de la condition physique.

En se basant sur l'apprentissage mé canique et les techniques d'apprentissage statistique, jusqu'à pré sent, plusieurs modè les de substitution ont é té utilisé s dans le calcul é volutionnaire assisté . Les modè les de substitution les plus populaires, y compris la ré gression polynomiale, le modè le de krigeage, les fonctions de base radiales, les ré seaux de neurones et les machines de vecteurs de soutien.

Modèle de krigeage

Le modè le de krigeage est une mé thode de modé lisation basé e sur le processus de Gauss pour interpoler des donné es dé terministes sans bruit et s'est ré vé lé utile dans une grande varié té de domaines. La caracté ristique distinctive du modè le de Krigeage est qu'il fournit non seulement une ré ponse pré dite (moyenne de pré diction) à un point non é chantillonné , mais aussi une estimation de la variance de pré diction (ou é cart type de pré diction). Cette variance donne une indication de l'incertitude dans le modè le de krigeage, qui ré sulte de la construction de la fonction de covariance. Cette derniè re est basé e sur l'idé e que lorsque les points d'entré e sont proches l'un de l'autre, la corré lation entre leurs sorties correspondantes sera é levé e. Par consé quent, l'incertitude associé e aux pré dictions du modè le sera faible pour les points d'entré e qui sont prè s des points de formation, et augmentera à mesure que l'on s'é loigne des points de formation. L'é cart-type de la pré diction et la moyenne de la pré diction fournissent des informations pré cieuses pour é quilibrer l'exploitation et donc bé né fique pour le calcul é volutionnaire assisté . Compte tenu de cet avantage, le modè le de krigeage est utilisé comme substitut dans la straté gie d'é volution assisté e par substitution. 3) KA-CMA-ES en utilisant la procé dure de classement approximatif (ARP). Compte tenu à le taux de ré ussite et le performance de l'accé lé ration, on peut conclure que l'ARP-EI est pré fé rable parmi tous les algorithmes KA-CMA-ES é tudié s dans cette thè se.

Quatrième

Figure 1 .

 1 Figure 1.1 Illustration of local and global minimum.......................................................................... 7

1 7 Figure 1 . 2 20 Figure 1 . 3 25 Figure 1 . 4 26 Figure 2 . 1 35 Figure 3 . 1 47 Figure 3 . 2 A 3 - 49 Figure 3 . 3 51 Figure 3 . 4 61 Figure 3 . 5 63 Figure 3 .

 1712201325142621353147323493351346135633 Figure 1.1 Illustration of local and global minimum.......................................................................... 7 Figure 1.2 Working styles of fitness approximation ........................................................................ 20 Figure 1.3 Illustration of fixed individual-based control.................................................................. 25 Figure 1.4 Illustration of fixed generation-based control ................................................................. 26 Figure 2.1 Gaussian mutation operators in evolution strategies....................................................... 35 Figure 3.1 Process of surrogate modeling ........................................................................................ 47 Figure 3.2 A 3-level and 3-dimensional full factorial design (27 points) ........................................ 49 Figure 3.3 LHS realization of 10 samples in a two-dimensional design space ................................ 51 Figure 3.4 Structure of RBF network models .................................................................................. 61 Figure 3.5 Structure of a three-layer MLP network ......................................................................... 63 Figure 3.6 The -insensitive loss function for SVR ......................................................................... 65

Figure 4 . 1

 41 Figure 4.1 Illustration of standard (,)-ES and the surrogate-assisted ES using pre-selection. ..... 81

Figure 4 . 2

 42 Figure 4.2 Individual-based generation-based evolution control in surrogate-assisted ES .............. 83

Figure 4 . 3 Figure 4 . 4 Figure 4 . 5 Figure 4 . 6 Figure 4 . 7 Figure 4 . 8 Figure 4 . 9 Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 4 .Figure 5 . 1 Figure 5 . 2 Figure 5 . 3 Figure 5 . 4 Figure 5 . 5 Figure 5 . 7 Figure 5 . 8 Figure 5 . 9 Figure 5 .Figure 5 .Figure 6 . 2 Figure 6 . 3 Figure 6 . 4 Figure 6 . 5 Figure 6 . 8 Figure 6 . 9

 43444546474849444444444444444515253545557585955626364656869 Figure 4.3 Radar charts of average success rate and speedup performance of KA-CMA-ES using preselection without model impact control (PS). ....................................................................... 116 Figure 4.4 Average success rate and outperform rate of pre-selection with model impact control using different model quality measures. ......................................................................................... 119 Figure 4.5 Average success rate and speedup of pre-selection without and with model impact control. ............................................................................................................................................... 122 Figure 4.6 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using fixed individual-based control with different metrics. .................................................................... 125 Figure 4.7 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using mixed individual-based control (MIC). ............................................................................................ 127 Figure 4.8 Average SR and SPU of KA-CMA-ES using fixed and mixed individual-based control (FIC and MIC). ...................................................................................................................... 130 Figure 4.9 Average success rate and speedup performance of KA-CMA-ES using ARP. ............ 131 Figure 4.10 Average success rate and speedup of KA-CMA-ES using FGC and AGC. ............... 135 Figure 4.11 Comparison of average success rate (SR) of KA-CMA-ES algorithms. .................... 138 Figure 4.12 Comparison of average speedup performance (SPU) of KA-CMA-ES algorithms. .. 139 Figure 4.13 Convergence graphs of function 1 .............................................................................. 140 Figure 4.14 Convergence graphs of function 2 .............................................................................. 141 Figure 4.15 Convergence graphs of function 3 .............................................................................. 141

1

 1 Critical values of model quality measures ....................................................................... 92Table 4.2 Test functions for experimental studies ......................................................................... 111 Table 4.3 Results of KA-CMA-ES using pre-selection without model impact control (PS) ......... 115 Table 4.4 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using preselection without model impact control (PS). ....................................................................... 116 Table 4.5 Results of KA-CMA-ES using pre-selection with model impact control (CPS). .......... 118 Table 4.6 Average success rate (SR) and speedup (SPU) performance of KA-CMA-ES using preselection with model impact control (CPS). .......................................................................... 119 Table 4.7 Comparison of KA-CMA-ES using pre-selection without and with model impact control, including PS-Interval, CPS-w Q , selection CPS-Q and rank CPS- . ................................................... 121 Table 4.8 Average success rate and speedup performance of KA-CMA-ES using pre-selection without and with model impact control, including PS-Interval, CPS- . .............................................................................................................................. 122 Table 4.9 Results of KA-CMA-ES using Fixed Individual-based Control (FIC) with different metrics. ............................................................................................................................................... 124 Table 4.10 Average success rate and speedup performance of KA-CMA-ES using fixed individualbased control (FIC) with different metrics. ........................................................................... 125

Algorithm 4 . 1

 41 Approximate Ranking Procedure ............................................................................. Algorithm 4.2 Pre-Selection Procedure............................................................................................ Algorithm 4.3 Fixed Individual-based Control (FIC) using Metric ................................................. Algorithm 4.4 Mixed Individual-based Control (MIC) .................................................................... Algorithm 4.5 Modified Approximate Ranking Procedure (ARP) .................................................. Algorithm 4.6 The KA-CMA-ES using Pre-Selection ................................................................... Algorithm 4.7 The KA-CMA-ES using Individual-based Control ................................................ Algorithm 4.8 The KA-CMA-ES using Approximate Ranking Procedure (ARP) ........................ Algorithm 4.9 The KA-CMA-ES using Adaptive Generation-based Control (AGC) ................... Algorithme 6.1 Le   , -CMA-ES w    ......................................................................................

  of d design variables, vectors l and u are called lower and upper bounds, respectively. Box constraints restrict the search space to the hyperrectangle  black-box property and complex landscape of objective function in expensive optimization problems, evolutionary algorithms (EAs), which is a class of derivative-free and powerful global optimizers, are appropriate for solving expensive optimization problems. However, the main difficulty in applying EAs to solve expensive optimization problems is that EAs usually need a large number of fitness function (objective function) evaluations before obtaining a satisfying result. Consequently, surrogate-assisted evolution algorithms were motivated from reducing computational costs in evolutionary optimization of expensive problems. Optimization using surrogate-assisted evolutionary algorithms is also known as surrogate-assisted evolutionary optimization or surrogateassisted evolutionary computation in the optimization community. This work concentrates on surrogate-assisted evolution strategy (ES) for expensive optimization problems. The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) and Kriging model are chosen as the two components of surrogate-assisted evolution strategy. This Kriging-Assisted Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is abbreviated as KA-CMA-ES in the thesis. Our goal is to investigate the existed surrogate-assisted ES and develop new efficient algorithms of KA-CMA-ES for expensive problems. A comprehensive study on KA-CMA-ES are performed in this thesis. New training set selection methods, pre-selection strategy, evolution control and approximate ranking procedure have been modified or developed, and the corresponding KA-CMA-ES algorithms are formulated. Then, experimental studies of KA-CMA-ES algorithms using existed training set select methods, pre-selection and evolution control, and our new developed KA-CMA-ES algorithms are carried out to analyze and evaluate the performance of the algorithms.

Thesis Outline 4 CMA

 4 -ES, which is one of the most successful evolution strategies and is chosen as the representative of ESs, is completely introduced. Chapter 3 introduces the surrogate modeling of computer experiments. Design of experiments, surrogate models (including polynomial regression model, Kriging model, radial basis function model, multi-layer perceptron networks and support vector regression), and model validation are successively described in this chapter. Chapter 4 dedicates to the Kriging-Assisted CMA-ES (KA-CMA-ES) algorithms, which combine the Kriging model and CMA-ES to reduce the computational cost of ES for expensive optimization problems. A brief survey of surrogate-assisted evolution strategies and some open issues in this topic are firstly provided. We focus on dealing with these open questions and developing new KA-CMA-ES algorithms for expensive optimization problems. Then, a series of KA-CMA-ES algorithms using pre-selection strategy, individual-based control, approximate ranking procedure and generation-based control are illustrated in mathematical formulations. Finally, experimental studies on KA-CMA-ES algorithms are performed on forty test problems (including twelve benchmarking functions),

1. 1

 1 Basic Concepts of Optimization 6 can be taken values. In other words, the values that design variables taken should satisfy the constraints, if constraints exist. In this work, only the continuous optimization problems with box constraints are considered. The standard mathematical form of an (continuous) optimization problem is called equality constraints. Conventionally, the standard form defines a minimization problem. A maximization problem can be treated by negating the objective function. The constraints on design variables define the set or domain

Figure 1 . 1

 11 Figure 1.1 Illustration of local and global minimum

1. 2

 2 Overview of Optimization Algorithms 8 classified as derivative-based and derivative-free, deterministic and stochastic, trajectorybased and population-based, and so on. In this section, the derivative-based and derivativefree algorithms classification is considered. Derivative-based methods (or gradient-based algorithms) use the derivative information in the search process, while the derivative-free methods (gradient-free algorithms) only use use the values of the objective, not any derivatives.Derivative-based methods generally require a much smaller number of iterations to converge to an optimum compared to derivative-free methods. However, only convergence to local minimum is guaranteed for derivative-based methods, while derivative-free methods are able to find global minimum. In addition, derivative-based methods are limited in cases where the objective function is always differentiable and has continuous derivatives over the search domain. Derivative-free methods do not place limitations on objective functions.

Algorithm 1 . 5

 15 Trust-Region Algorithm 1: input: objective function   f x , initial point   0 x , initial trust-region radius   0

Algorithm 1 . 6 7 :

 167 General Framework of an Evolutionary Algorithm until termination criterion fulfilled

Figure 1 . 2

 12 Figure 1.2 Working styles of fitness approximation

,

     new offspring individuals are firstly created from  parents. Then,  individuals out of these  individuals are pre-selected according the merit function values evaluated by Gaussian process model to generate the offspring population. The pre-selected  offspring are evaluated by original fitness function for selection and recombination.

  are two approaches to fixed evolution control: individual-based evolution control and generation-based evolution control. Individual-based evolution control means that in each generation, some of the individuals are evaluated by surrogate model and the others are evaluated using the original fitness function, which is illustrated in Figure 1.3. In individualbased evolution control, the individual selection can be random or using some strategy, e.g., selecting the best individual (based on the prediction of surrogate model) for evolution control. In generation-based evolution control, all individuals in a selected generation will be evaluated by the original fitness function. The generation selection can be random or with a fixed frequency. A generation-based evolution control is shown in Figure 1.4.

Figure 1 . 3 26 Figure 1 . 4

 132614 Figure 1.3 Illustration of fixed individual-based control



  evolution strategies, the commonly used selection strategies are the comma (indicated by ,) selection and plus (indicated by +) selection:  Plus selection selects the  -best solutions from the union of the last parent population (  individuals in parent population) and the current offspring population (  individuals in offspring population), and is correspondingly Comma selection, i.e.,   , -ES  , selects -best  solutions exclusively from the offspring population (  individuals in offspring population), neglecting the parent population, even if the parents have a superior fitness.

Figure 2 . 9 )

 29 Figure 2.1 (a).

  distribution. For a problem of d-dimensional search space, the number of endogenous strategy parameter in case of Equation (2.8) is   1 O , i.e., constant. In case of Equation (2.9), a vector of endogenous parameters with size   Od is required. To adapt an arbitrary covariance matrix C in Equation (2.10),   2 Od endogenous parameters are required [3].

Figure 2 . 1

 21 Figure 2.1 Gaussian mutation operators in evolution strategies

  with Cumulative Step-Size Adaptation (CSA).



  denotes the step-size (overall standard deviation) at generation g . Its initial values   0

  -size  and c d  p for covariance matrix C adaptation, are used in CMA-ES.

1 c

 1 and c  are the learning rate for rank-one-update and rank--update, respectively. Their default values are:

Section 3 . 1 ,

 31 the general process of surrogate modeling is introduced. The three steps of surrogate modeling, (i) design of experiments; (ii) training of the surrogate model; and (iii) model validation, are respectively addressed in Section 3.2, 3.3 and 3.4.

  represents the input-output relationship exhibited by simulation model over a given design domain." The set of sampling points X is determined by design of experiments (DOE), and the corresponding responses y are generated by the experiments conducted under that DOE. Thus, design of experiments is the first step in the construction of surrogate models. The samples (sample points) used to construct the surrogate models are generally called training points. The set of data that includes input data and the corresponding output data is called training data set, while the construction of the surrogate models based on training data set is often called model training. With the training data set, a variety of surrogate models can be trained. Before the trained surrogate model is used in lieu of the original simulation model, it is necessary to evaluate the performance or expected accuracy of the surrogate model. Therefore, the typical process of surrogate modeling mainly involves three steps (i) design of experiments; (ii) construction or training of the surrogate model; and (iii) model validation. This process is illustrated in Figure 3.1. These three issues for surrogate modeling are discussed separately in following sections.

Figure 3 . 1

 31 Figure 3.1 Process of surrogate modeling

Figure 3 .

 3 Figure 3.2 A 3-level and 3-dimensional full factorial design (27 points)

  Designs generated by Latin hypercube sampling are called Latin hypercube designs (LHD) in computer experiments design. Latin hypercube sampling (LHS) was proposed by McKay et al. (1979) specifically for computer experiments in which inputs are chosen randomly from some specified distribution. Latin hypercube designs have one-dimensional uniformity in that, when projected on each dimension, each interval of the dimension has exactly one observation.

Figure 3 . 3

 33 Figure 3.3 LHS realization of 10 samples in a two-dimensional design space

3 .

 3 Surrogate Modeling of Computer Experiments 57 Depending on the form of the regression function, Kriging has been prefixed with different names. Simple Kriging assumes the regression function to be a known constant. A more popular version is ordinary Kriging, which assumes a constant but unknown regression function. In general, universal Kriging treats the trend function as a polynomial regression model, namely,

Zx,

  between any two points x and  x , xx is the correlation function between x and  x , and 2  is the stochastic process variance. The correlation function in Equation (3.16) affects the both the range of influence and the smoothness of the model [55]. The most commonly used correlation function is Gaussian correlation function, which is defined by unknown parameters of the correlation function and can be contained in a column vector d is the dimension of design space (number of design variables), k x and k x denote the k-th component of the samples x and  x , respectively, and kk xx   is the absolute distance between k x and k x . Consider on the training data set

  From above description, in order to build Kriging surrogate model, the regression coefficients β in Equation (3.15), the correlation parameter θ in Equation (3.17), and the stochastic process variance 2  in Equation (3.16) need to be determined. The maximum likelihood estimate method is used to determine the unknown model parameters β , 2  and θ .

  (3.25) and(3.26). So far, we have determined the values of all unknowns in Kriging model, i.e., the Kriging model has been completely built. The prediction at any untried point can be estimated by Equation(3.18).

w

  where x is the input, m is the number of neurons in hidden layer,   i c is the center or prototype of the i-th radial basis function   i  x , and  is usually a Euclidean norm. For input the output of the RBF network   f x is given by a linear combination of a set of radial basis functions is the weight for i-th radial basis function. The RBFs im  transform the input in terms of the distance between the input x and the center   i c . For the generalized RBF network, the center   i c are also unknown and have to be learned by other methods such as the k-means method[START_REF] Tenne | Computational Intelligence in Expensive Optimization Problems[END_REF].



  Inverse multi-quadric:     12 22 rr    .Whether we choose a set of parametric basis function or non-parametric ones, the weights via the interpolation condition[START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF]. With the training points

Figure 3 . 4

 34 Figure 3.4 Structure of RBF network models

Figure 3 . 5

 35 Figure 3.5 Structure of a three-layer MLP network

  Note that the tacit assumption in Equation(3.35) is that a function   f x exists that approximates all pairs set with precision  , or in other words, that the convex optimization problem in(3.35

3 .Figure 3 .

 33 Figure 3.6 The -insensitive loss function for SVR

42 )

 42 Substituting Equations (3.39), (3.40), (3.41) and (3.42) into Equation (3.38) to eliminate   i  and   , i   finally yields the dual optimization problem:

  model and that predicted by the surrogate model. Let of input variables, y denotes the output from simulation model   f x , and ŷ signifies the output predicted by surrogate model   f x built on training data set

3. 4

 4 Model Validation/Model Performance Assessment 68

2 R

 2 m is the number of training pints (sampling points),   , i y   ˆi y and y respectively represent the (or experimental) response, the prediction and the mean of responses. The value of coefficient of determination 2R is between 0 and 1, i.e., is near to 1, generally, it indicates the model fit well the sampled data.

  optimization procedure based on the combination of ES and ANN was found to be very effective in shape and sizing structural optimization problems.Individual-based and generation-based evolution control methods were proposed in[START_REF] Jin | Managing approximate models in evolutionary aerodynamic design optimization[END_REF] for surrogate-assisted evolution strategies with CMA-ES. Evolution control is popular in surrogate-assisted evolutionary algorithms, including evolution strategies. In individualbased control, part of the individuals in the population are chosen and evaluated with the original fitness function. The controlled individuals can be chosen by random or best strategies. In generation-based control, the whole population of  generations will be evaluated with original fitness function in every p generations, where p   . Reference[START_REF] Jin | Managing approximate models in evolutionary aerodynamic design optimization[END_REF] focused on generation-based evolution control and proposed a strategy to determine the control frequency p  . Specifically, the current model error was used to estimate the local fidelity of the approximate model and the local model fidelity was adopted to determine the control frequency. After each generation control, new data points were added to training data set and then the approximate model was updated. In the on-line learning of approximate model, the weighted learning using covariance matrix was proposed and used in this work.

  Neighbor Points Selecting evaluated points that are more relevant to current search space is reasonable. Considering that ES use Gaussian mutation, the k-nearest neighbor points ( training kn  ) to the mean of the Gaussian distribution m best represent the actual search space and are good candidates for model training. The Euclidean distance is usually used as similarity measure metric. However, when f is a multi-modal function, it may happen that m is located around a local optimum which has been visited by the evolutionary search many generations ago, and many of training n closest points to m are correspond to a smaller cluster of this local optimum.

  correlates with the ability to correctly select individuals in evolutionary computation, and selection based or ranking based model quality measurements were suggested to use. From the perspective of evolutionary computation, the correct selection is of importance. In [82], Jin et al. proposed model quality measures that focus primarily on the correct model-based selection rather than the approximation accuracy of surrogate model. Some of these measures are presented below. Let us consider the case of the (,)-selection. The model-based selection process (selection based on surrogate model) selects  out of the  individuals with the best predicted fitness. An individual is correctly selected, if it would also be selected by a selection process based on the original (true) fitness of the individuals. The first selection based measure is about the number of individuals that have been selected correctly using the surrogate model: select ,

  et al.[START_REF] Jin | Quality measures for approximate models in evolutionary computation[END_REF] suggested to use the rank correlation coefficient[START_REF] Spearman | The Proof and Measurement of Association between Two Things[END_REF] (also known as Spearman's rank correlation coefficient), given by of the monotonic relation between the ranks of two variables, to evaluate the quality of the surrogate model. In the case of using rank  to measure the monotonic relation between the individual's rank based on surrogate model and that based on original fitness function, i d is the difference between the ranks of the i-th offspring (among  ) individual based on the original fitness function and on the surrogate model. The rank  has the range rank 11     . The higher the value of rank  , the stronger the monotonic relation with a positive slope between the ranks of the two variables. When the model predicts correct ranking of all  individuals, we have rank 1   . Inversely, when the model predicts inverse ranking of individuals, rank 1   . The rank  calculated by Equation (4.6) takes the ranks of all individuals rather than only take the ranking of the selected individuals into account.

2 . 1 .

 21 As previously description, in pre-selection strategy, for each generation, Pre   individuals are generated (through mutation) and evaluated by surrogate model f , then  out of Pre  best individuals are selected (based on the fitness value estimated by f ) to re-evaluate by the original fitness function f . Thus, all the  parent for next generation selected from  offspring are evaluated by real fitness function f . The pre- selection strategy can be illustrated in Figure 4.1. Evolution control takes place when surrogate model f is used in fitness evaluation of ES. In individual-based evolution control,  individuals are firstly evaluated using the surrogate f , then    individuals are selected to re-evaluated using the original fitness function f . As a result, the fitness values of  parent for next generation can be based on surrogate f or on real fitness f . In generation- based evolution control, evolution control is carried out at a specific frequency, where the whole population of the controlled generation are evaluated using the original fitness function f , while in other generations the whole population are evaluated by surrogate ˆ. f Surrogate-assisted ES using individual-based and generation-based evolution control are illustrated in Figure 4.2. Some open issues in pre-selection and evolution control are described in the following.

Figure 4 . 1

 41 Figure 4.1 Illustration of standard (,)-ES and the surrogate-assisted ES using pre-selection.



  was dynamically controlled by model quality measurement based on the number of correctly pre-selected individuals. It was demonstrated in[START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF] that controlling the impact of model impact ( Pre  ) enhances the performance of model assisted Evolution Strategies (MAES) with fixed model influence (fixed Pre  ). The concept of model impact control will attract much more attention in MAES. Another issue in pre-selection is the pre-selection criterion which is used to identify the most promising individuals. For Gaussian process (Kriging) model assisted ES, several preselection criteria have been applied in literature. For other surrogate model assisted ES, usually, the model prediction is used to pre-screen the most promising individuals. However, no precise conclusion on performance of different pre-selection criteria has been drawn so far. Comprehensive study of existed pre-selection criteria is still needed. And more sophisticated pre-selection criterion may be developed in the future. In Evolution Control, the main issues is to determine the control frequency, which denotes the number of individuals controlled in individual-based control and frequency of controlled generation in generation-based control. The control frequency plays important role in guaranteeing the correct convergence of the surrogate-assisted evolutionary optimization. Specifically, for individual-based evolution control, in each generation  individuals are selected from  (  is the size of population of the standard ES) offspring for control purpose, i.e, to be evaluated by original fitness function. If the controlled individuals are chosen randomly, it is called a random strategy. If the best  individuals are selected to be controlled, we called it a best strategy. In generation-based evolution control, the whole population of  generations will be evaluated with original fitness function in every p generations, where p   . Then, the control frequency in individual and generation-based evolution control can be presented by the fraction and p

4 .

 4 Kriging-Assisted CMA Evolution Strategy 83[START_REF] Jin | Managing approximate models in evolutionary aerodynamic design optimization[END_REF] introduced a framework to adapt the control frequency based on the fidelity of the surrogate model.

Figure 4 . 2

 42 Figure 4.2 Individual-based generation-based evolution control in surrogate-assisted ES



  the actual (current) search subspace (the region that the distributed offspring can reach) is determined by current distribution of offspring   , d   . It is apparent that the sampled offspring probably lie in the region with high probability of the multivariate normal distribution. Consequently, the region yielded by the confidence interval of the multivariate normal distribution can be used to represent the actual subspace with a given probability. The confidence interval for the multivariate normal distribution is the quantile function for probability p   01 p  of the chi-squared distribution with d degrees of freedom. The inequality (4.8) provides the confidence region containing p of the probability mass of the MVN distribution   , d   . In other words, a random vector from distribution   , d   has probability p of satisfying the inequality (4region defined by inequality (4.8) has a probability p that it cover the offspring. Put another way, this region is the subspace that the offspring generated from the MVN distribution   , d   can reach with a probability p . The probability p indicates the extent that the confidence region covers the actual search subspace. Therefore, we propose to select the training data set training by

4. 3 . 2 4 . 3 . 2 . 1

 324321 Pre-Selection Strategy with Model Impact Control Introduction of Pre-Selection with Model Impact Control (CPS)The size of pre-selection population size Pre  controls the impact of the surrogate model f on the evolutionary optimization process. Thus, the search procedure would benefit from the appropriate control of Pre  . Simply, the pre-selection population size is set as a constant by users during the optimization. Recently, researchers attempt to adapt Pre  during the evolutionary search process. It is intuitive to use the model quality to control the preselection population size Pre  . In[START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF], the value of Pre  was controlled by model quality measurement based on the number of correctly pre-selected individuals. This idea gives significant guidance for us in pre-selection population size control.In this paragraph, we briefly present the model quality measurement and adaptation of Pre that were proposed by Ulmer et al. [84]. Since in optimization by evolution strategies, only the correct selection is of importance, the authors used a selection based model quality which was similar to that proposed by Jin et al. [82]. Consider the model based selection process selects  out of  individuals with the best predicted fitness. An individual is 4. Kriging-Assisted CMA Evolution Strategy 87 correctly selected, if it would also be selected by a selection process based on the original (true) fitness of the individuals. Thus, the number of correctly selected individuals can be used to measure the quality of the model based selection process. The authors assumed that the selection quality determined for a selection of  out of  individuals, for which we known the predicted and the true fitness, is equivalent to the selection quality of the pre-selection. If the individual with the i-th best true fitness is correctly selected, a rank (grade) of   i   is given to this individual (this rank could be considered as the grade of the individual in selection). The summed rank of all correctly selected individuals is defined as the quality selection Q of the model based selection process. The difference between selection Q and  in Equation(4.4) is that selection Q only counts the grades of correctly selected individuals but counts the grades of all selected individuals including both correctly and incorrectly selected individuals. When none individual is correctly selected, selection Q has its minimum value all individuals are correctly selected, the maximum model quality is

  i.e., dynamically controlling Pre  based on model quality. In this thesis, a variant of the selection based model quality measurement, which was proposed by Ulmer et al. in[START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF] and described above, is defined based on the weights used in recombination operator in an evolution strategy.We consider an evolution strategy that use weighted recombination operator, recombination, the weight values depend on the fitness ranking, in that better parents get larger weights than inferior ones. The widely used superlinear decrease weights are given by

.

  recombination weights calculated by real fitness ranking is assigned to this individual as the score it obtains in model selection; otherwise, it has score zero. The measurement of selection based model quality w Q is defined as the summed score of all model selected individuals. Obviously, a larger value of w Q means the surrogate model has higher quality of correctly selecting individuals, and vice versa. If the model selects none correct individual, In other cases, i.e., only part of model selected individuals are correct, 01 w Q  . From the definition of w Q , it can be found that when 4.3 Address Some Open Issues of Surrogate-Assisted ES 90 individuals with higher real fitness ranking and/or more individuals are selected correctly, w Q obtains higher value and the model quality is better. Similar to Equation (4.11), the expectation value of w Q for a purely random selection process is simply given as the product of the expectation value of the number of correctly selected individuals and the expectation value of the score of a correctly selected individual :

2 . 1 .C

 21 In pre-selection procedure, Pre   individuals are generated, then all Pre  individuals are evaluated by surrogate model f and the estimated fitness values are used to pre-select the  best individuals, which will be evaluated by the original fitness function. After the initial sampling, the pre-selection procedure can be performed in each iteration of the evolution loop. The module of pre-selection for evolution loop is expressed in Algorithm 4.2. for re-evaluation using original fitness function In Algorithm 4.2, the current distribution mean vector   g m , step-size   g ,and the archive which contains all previously evaluated data points are used in pre-selection (Line 1). The pre-selection population size for current generation   Pre g  is the parameter of pre-selection procedure. With the inputs, training set T is firstly selected by the training set selection method (Line 2), which can be the recently evaluated points, knearest neighbor points to distribution mean and confidence interval method proposed in Section 4.3.1, and the surrogate model f is trained from the selected training set (Line 3).

4. 3

 3 Address Some Open Issues of Surrogate-Assisted ES 94 function is called a controlled generation. Correspondingly, evolution control can be divided into two approaches: individual-based and generation-based evolution control. In individualbased evolution control, a certain number of individuals within a generation are selected to be evaluated using the original fitness function, and other individuals in the generation use the surrogate model for fitness evaluation. Generation-based evolution control is carried out at a certain frequency. In a controlled generation, all individuals are evaluated with the original fitness function. The most important question for evolution control, both of individual-based and generation-based evolution control, is how to determine the control frequency, which signifies the number of controlled individuals and generations for individual-based and generation-based evolution control, respectively, in order to guarantee the correct convergence of the evolutionary search when false optima are present in the surrogate model of fitness function. For individual-based evolution control, both fixed and adaptive control methods can be used. In fixed individual-based evolution control, the number of controlled individuals is fixed, i.e., a fixed number of individuals are selected to re-evaluated by using the original fitness function. While, for adaptive individual-based evolution control, the number of reevaluated individuals is adaptive, usually depend on the fidelity (quality) of the surrogate model. The controlled individuals can be selected randomly or by a best strategy, which selects best individuals based on predicted fitness. Obviously, fixed individual-based control is simpler and thus more easily to implement than adaptive control. However, adaptive individual-based evolution control may require less number of re-evaluated individuals

x 4 :n

 4 select the init n best individuals based on metric  selected individuals by   f x and add to the set , of selecting only one individual for re-evaluation and updating the surrogate before entering the iteration loop of approximate ranking procedure, we select individuals to be re-evaluated by   f x for updating the model. Obviously, more information would be added through this change. Secondly, in the iteration loop of approximate ranking procedure, we use a batch size

  where

  vectors (search points, candidates or individuals), and dd   C is the current covariance matrix of the CMA-ES. This covariance matrix based distance takes into account the correlations of the data. In this work, a transformation of variables which transforms correlated variables into uncorrelated is suggested to be performed before the model training and prediction. The CMA-ES, which uses general (or correlated) mutation operator, is considered. In each generation of the CMA-ES,  offspring are generated by mutation

  Therefore, each time before training a surrogate model and making predictions by the model, we firstly transform all the input vectors of previously evaluated points and the current  offspring individuals using Equation (4.24). Then, the training set is selected by using the k-Nearest Neighbor Points based on Mahalanobis Distance, the confidence interval of MVN or other methods. That is to say, model learning and prediction are performed on transformed variables.

4. 4

 4 Kriging-Assisted CMA-ES (KA-CMA-ES) Algorithms 104In other words, before the evolution loop of KA-CMA-ES, an initial design of experiments (DOE) is carried out by using Latin Hypercube Designs (LHD) based on maximin distance criterion. As for the size of the DOE, DOE n , the DOE 10 nD  rule of thumb which was studied in[START_REF] Loeppky | Choosing the Sample Size of a Computer Experiment: A Practical Guide[END_REF] is adopted. Specifically, DOE 10 nD  ( D is the dimension of the problem) sampling points are generated by LHD. However, in KA-CMA-ES using pre-selection, the size of preselection population is set as Pre 2   . In order to runs pre-selection after the initial sampling, DOE n should not smaller than Pre  . Thus, we correspondingly set CMA-ES algorithms which are expressed and investigated subsequently.

Algorithm 4 . 6

 46 without model impact control (PS) and pre-selection with model impact control (CPS) are incorporated into CMA-ES. The standard   , -CMA-ES w    , which has already presented in Section 2.3.3, is not repeated here. The KA-CMA-ES using pre-selection can be easily developed by introducing initial sampling and informed start point before the evolution loop, and then runs the pre-selection procedure before the evaluation of offspring population (using original fitness function f ) in the standard evolution loop of the   The KA-CMA-ES using Pre-Selection 1: given: strategy parameters of CMA-ES, pre-selection method (PS or CPS) and parameters ( Pre  model quality measure Q for CPS).

  as standard CMA-ES algorithm Line 9~13) 24: 1 gg  25: until termination criterion is fulfilled In each evolution iteration, the pre-selection procedure is firstly performed. After the pre-selection step, the evaluation, selection and recombination are performed on the  pre- screened offspring in the same way as in standard CMA-ES. If there is no model impact control in pre-selection (PS), the size of pre-selection population Pre  is set by user and keeps constant during the computation. For pre-selection with model impact control (CPS), the model quality Q (including w Q , selection Q and rank  ) need to be estimated after the evaluation of  offspring by original fitness function, and then the size of pre-selection population Pre  is updated according to Equation (4.16). The pseudocode of KA-CMA-ES using Pre-Selection without and with model impact control is given in Algorithm 4.6. This algorithm contains two methods: KA-CMA-ES using PS, which presents Kriging-assisted CMA-ES using preselection without model impact control, and KA-CMA-ES using CPS,which stands for Kriging-assisted CMA-ES using pre-selection with model impact control. The pre-selection procedure used in Line 10 and 12 has been given in Algorithm 4.2.

ForAlgorithm 4 . 7

 47 Kriging-Assisted CMA-ES using individual-based evolution control, fixed individual-based control (FIC) using metric and our proposed mixed individual-based control (MIC) can be used. In fixed individual-based control (FIC), metric is used for selecting the most promising individuals for re-evaluation. The details about metrics (Mean, SD, SLB, POI and EI) for FIC have been discussed in Section 4.3.3.1. In mixed individualbased control (MIC), two metrics, i.e., the Mean (mean of prediction) and SLB (the statistical lower bound), are provided. The pseudocode of KA-CMA-ES using individual-based control, including FIC and MIC, is presented in Algorithm 4.7. The initial sampling and informed start point are firstly introduced in KA-CMA-ES using individual-based control (Line 3~5), before the evolution loop. In each evolution iteration,  individuals (offspring) are generated by mutation operator (Line 7~10). Then, the fixed individual-based control (FIC) or mixed individual-based control (MIC) is performed 4. Kriging-Assisted CMA Evolution Strategy 107 (Line 13~17). With the results of FIC or MIC, the distribution mean, step-size and covariance matrix can be updated (Line 18) the same as the standard CMA-ES. The KA-CMA-ES using Individual-based Control 1: given: strategy parameters of CMA-ES, individual-based control method (FIC or MIC) and parameters (  for FIC, min max , for MIC), and metric   C x used in individual control.

20 :

 20 until termination criterion is fulfilled 4.4.5 KA-CMA-ES using Approximate Ranking Procedure (ARP) Previously described modified approximate ranking procedure, which has been illustrated in Algorithm 4.5, is embedded into CMA-ES. The pseudocode of KA-CMA-ES 4.4 Kriging-Assisted CMA-ES (KA-CMA-ES) Algorithms 108 using approximate ranking procedure (ARP) is given in Algorithm 4.8. The initial sampling and informed start point are firstly performed (Line 3~5). In each evolution loop, the approximate ranking procedure is called after  individuals have been generated.

Algorithm 4 . 8

 48 The KA-CMA-ES using Approximate Ranking Procedure (ARP) 1: given: strategy parameters of CMA-ES, parameters for approximate ranking procedure ( init b

14 :Algorithm 4 . 9

 1449 until termination criterion is fulfilled4.4.6 KA-CMA-ES using Adaptive Generation-based Control (AGC)The CMA-ES is further combined with the adaptive generation-based control (AGC) strategy described in Section 4.3.5, in which the model quality is used to control the frequency of generation control. Specifically, in adaptive generation-based control, if the model quality   g Q of current controlled generation exceeds the threshold value TS Q , next generation is evaluated using the surrogate model; otherwise, next generation is controlled. The pseudocode of KA-CMA-ES using adaptive generation-based control (AGC) is expressed in Algorithm 4.9. The model quality measures w Q , selection Q and rank  can be used in AGC. In Algorithm 4.9,   control g I indicates the generation g is controlled or not. If   control 1 g I  , generation g is controlled (evaluated by the original fitness function), otherwise, generation g is evaluated by the surrogate model. The KA-CMA-ES using Adaptive Generation-based Control (AGC) 1: given: strategy parameters of CMA-ES, model quality Q for AGC 2

f

  x is objective or fitness function, D is the dimension of the problem (the number of parameters or variables), and S is the search space (search domain) defined by the lower bounds LB

4. 5

 5 Experimental Studies 112 CMA-ES, the initial mean   0 m of the evolution loop is the best candidate in initial sampling (informed start point), where by LHD within the search space. A single run of the algorithm is terminated, when the target function (fitness) one of the following termination conditions is satisfied: maxFES: maximum number of exact fitness function evaluations, for CMA-ES, stop if the range of the best objective function values of the last 10 30 D     generations and all function values of the recent generation is below 10 TolFun 10   . TolX: stop if all components of   g c p and all square roots of diagonal components of   g For each run of the problems, the number of exact function evaluations (FES) and current best fitness function values are recorded. All the best fitness function values should be from solutions that are evaluated using original fitness function. In other words, in KA-CMA-ES, the number of exact function evaluations (evaluating using original fitness function) and corresponding current best function values of evaluated individuals are recorded.

Furthermore, in order

  to analyze the performance of different training set selection methods, we have classified the total 40 test problems into different categories (or groups) according to their dimensionality and modality, and then evaluate the performance of PS using different training selection methods on each category and overall problems. The total test problems are grouped into: 2 dimensional problems (D=2), 5 dimensional problems (D=5), 10 dimensional problems (D=10), 20 dimensional problems (D=20), unimodal problem, and multimodal problems. For each category of problems and the overall problems, we compute the average success rate (SR) and average speedup performance (SPU). This is also used in subsequent investigation of other KA-CMA-ES algorithms. The average success rate and speedup performance are listed in Table 4.4 and plotted in radar charts (also known as spider chats or start plots) of Figure 4.3.

Figure 4 . 3

 43 Figure 4.3 Radar charts of average success rate and speedup performance of KA-CMA-ES using preselection without model impact control (PS).

From the experimental study

  of KA-CMA-ES using pre-selection without model impact control (PS), it can be found that the proposed Confidence Interval method for training set selection in KA-CMA-ES is appropriate. Considering the success rate and success performance as a whole, it is can be concluded that the 'Interval' (Confidence Interval) method for training set selection has better performance than 'Recently' and 'kNN' methods. By all accounts, the proposed Confidence Interval method clearly superior to the commonly used 'Recently' and 'kNN' methods for training set selection. Consequently, Confidence Interval method of training set selection is adopted for training set selecting in the remaining KA-CMA-ES algorithms owing to its superiority.4.5.2.2 Pre-Selection with Model Impact Control (CPS)In KA-CMA-ES using pre-selection with model impact control (CPS), three model impact control approaches, in which three different model quality measures including w Q , selection Q and rank  are used, are investigated. In experiments of KA-CMA-ES using CPS, we set the initial pre-selection population size  

wQQQ

  indicates pre-selection with model impact control using w Q as model quality measure, selection CPS-Q represents pre-selection with model impact control using selection Q as model quality measure, and rank CPS- signifies pre-selection with model impact control using rank  as model quality measure.According to the values of SP and SPU in Table4.5, for most problems, CPSperform worsen than CMA-ES on the problems of 20 dimensional 3 f and 4 f and 2 dimensional 10 f . Thus, it can be said that rank CPS- is more stable and reliable than CPSgenerally has better success performance and speedup performance.

4. 5

 5 Experimental Studies 118

Figure 4 . 4

 44 Figure 4.4 Average success rate and outperform rate of pre-selection with model impact control using different model quality measures.The average success rate and speedup performance of CPS are listed in Table4.6 and plotted in Figure4.4. For problems with D=2 and D=20, unimodal problems and multimodal problems, the

  for different categories of problems are also presented in Table4.[START_REF] Papadrakakis | Survey of quasi-Newton methods with reduced storage[END_REF] and plotted inFigure 4.5. In the radar chart of speedup performance in the right of Figure 4.5, the CPSw Q obviously performs better than other pre-selection methods. It is clear that the pre-selection with model impact control outperforms pre-selection without model impact control. All the pre-selection methods have average success rate larger than 0.9. The rank CPS- generally has the highest and stable success rate for all categories of test problems, especially on 20 dimensional problems and multimodal problems. Yet the average success performance (or speedup performance) of rank CPS- is not so good as that of CPSw Q . As discussed above, taking the success rate, speedup performance and the reliability into accounts, rank CPS- is preferable in all the investigated algorithms of KA-CMA-ES using pre-selection strategy. If only the success performance or speedup performance is considered, CPSw Q is outstanding in the investigated algorithms of pre-selection.

Figure 4 . 5 4 . 5 . 3 4 . 5 . 3 . 1 Fixed

 454534531 Figure 4.5 Average success rate and speedup of pre-selection without and with model impact control.

Figure 4 . 6

 46 Figure 4.6 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using fixed individual-based control with different metrics.

Figure 4 . 7 .

 47 In these tables, MIC-Mean represents mixed individual-based control using Mean as metric, and MIC-SLB means using SLB as metric in mixed individual-based control.

Figure 4 . 7

 47 Figure 4.7 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using mixed individual-based control (MIC).

Figure 4 .

 4 8, generally, KA-CMA-ES using FIC have higher success rate than that using MIC. The FIC-Mean and MIC-Mean has the lowest success rate among all the investigated algorithms using individual-based control. The FIC-EI generally has the highest success rate. Other algorithms have moderate average success rate larger than 0.85. However, KA-CMA-ES using MIC significantly outperform FIC according to the speedup performance and success performance. Particularly, MIC-SLB has the highest average SPU and acceptable success rate. Therefore, MIC-SLB is preferable among algorithms of KA-CMA-ES using individualbased control.

Figure 4 . 8

 48 Figure 4.8 Average SR and SPU of KA-CMA-ES using fixed and mixed individual-based control (FIC and MIC).

Figure 4 . 9

 49 Figure 4.9 Average success rate and speedup performance of KA-CMA-ES using ARP.

4. 5 . 6

 56 Performance Analysis of All the Investigated KA-CMA-ES Algorithms In previous three subsections, the KA-CMA-ES using pre-selection, individual-based control, approximate ranking procedure and generation-based control are separately investigated and analyzed. This subsection dedicates to comparing all the investigated algorithms of KA-CMA-ES. Based on previously discussion and results, only the best algorithms of each category of algorithms are chosen as the representatives and used in analysis and comparison here. The PS-Interval is taken as the representative of pre-selection without model impact control (PS). The CPSw Q and rank CPS- are considered as typical algorithms of pre-selection with model impact control (CPS). The FIC-EI is used to represent the fixed individual-based control (FIC). The ARP-EI stands for the KA-CMA-ES using approximate ranking procedure (ARP). For generation-based control, AGCare regarded as the representatives.

Figure 4 .

 4 Figure 4.11 and Figure 4.12.

Figure 4 .

 4 Figure 4.12 Comparison of average speedup performance (SPU) of KA-CMA-ES algorithms.

Figure 4 . 13 Convergence graphs of function 1 Figure 4 . 14 Convergence graphs of function 2 Figure 4 .

 413141424 Figure 4.13 Convergence graphs of function 1

Figure 4 . 16 Convergence graphs of function 4 Figure 4 .

 41644 Figure 4.16 Convergence graphs of function 4

Figure 4 .Figure 4 .

 44 Figure 4.18 Convergence graphs of function 6

Figure 4 . 20 Convergence graphs of function 8 Figure 4 .

 42084 Figure 4.20 Convergence graphs of function 8

Figure 4 . 22 Convergence graphs of function 10 Figure 4 .Figure 4 .

 4221044 Figure 4.22 Convergence graphs of function 10

D

  the overall area of the element defined by this normal vector n , and S is the effective resisting area. for metals) corresponds to the rupture of the element [102].

5 . 4 )

 54 Applications in Material Parameter Identification 155where  is the total strain tensor, the tensor e  and p  are known, respectively, as the elastic strain tensor and the plastic strain tensor. Correspondingly, the rate form of the additiveThe elastic law (coupled with damage) is expressed as: isotropic elasticity tensor,  is stress tensor and D is the scalar damage variable.Since von Mises yield criterion is used, the yield function for elastic-plastic damage model is given by q is the von Mises equivalent stress and y  is the yield stress which is defined by Swift hardening law as: and n are material parameters.The evolution of plastic strain is expressed by

17

 17

  )

  an elastic step and the elastic trial state coincides with the updated state at 1 n t  , i.e., we set     trial 11 nn     . Otherwise, we turn to the return-mapping step.

  To start with, let us consider the deviatoric/volumetric split of the elastic strain of Equation (5.18) 1 . This gives

D

   are unknown, other variables have defined previously. This nonlinear equation system is solved using Newton-Raphson method. There are two types of coupling can be used: strong coupling and weak coupling. The strong coupling solves the two equations simultaneously to obtain the two unknowns   and 1 n D  . While, the weak coupling firstly solves the first equation of the system (5.27) in which we assume that 1 n D  is equal to n D , in order to obtain   . Then the second equation is solved to compute the new damage value 1 n D  .

  strong coupling gives more accurate results, the weak coupling is generally preferred for computation time reasons in real simulation [103]. Hence, in this paper weak coupling is used to solve the return-mapping equation system (5.27). After obtaining   and 1 n D  , stress and other state variables can be updated. The complete stress 5. Applications in Material Parameter Identification 161 update algorithm for the numerical implementation of elastic-plastic damage model by means of the fully implicit elastic predictor/return-mapping scheme is summarized in Table

  Using numerical techniques to simulate the response of the structure with the given constitutive model and model parameters is called direct problem. While the parameter identification (or calibration), which determines model parameters on the basis of experimental data, is known as inverse problem [104]. The direct and inverse problems are sketched in Figure 5.1.

Figure 5 . 1

 51 Figure 5.1 Illustration of (a) direct problem and (b) inverse problem

Figure 5 . 2 .

 52 Figure 5.2. A general optimization problem for the inverse method of parameter identification can be expressed as

Figure 5 . 2

 52 Figure 5.2 Framework of inverse method for parameter identification

. 29 )p

 29 where p is the set of material parameters, are the numerical simulation responses with the input parameter of p corresponding to displacement D i , exp F i are the experimental responses corresponding to D i , and n is the number of data points. The second line of Equation (5.29) is the discrete formulation of the first line. In evaluation of the objective by its second line (discrete formulation), the simulation and experimental force-displacement data are interpolated at the same query points ( 12 0, D , D , , D n ) and then   f p is computed by the discrete formulation.

Figure 5 . 3

 53 Figure 5.3 Illustration of the difference between experimental and simulation results

  al.[108], this problem has been solved by adapting the integration domain of objective function in Equation(5.29). Let sim fracture be deal with. For both case, the objective function is evaluated between D0  and are completed with one new point defined as: the breaking force (load) remains the same whereas the displacement value is set as illustrated by the green line in the left of Figure 5.4. are completed with one new point defined as: the breaking force (load) remains the same whereas the displacement value illustrated by the green line in the right of Figure 5.4. With above adaptation of force-displacement data, the objective function (Equation (5.29)) can be obtained by evaluating the integration in Equation (5.29) using its discrete formulation.

Figure 5 . 4

 54 Figure 5.4 Tensile test force-displacement curve, Left: simulation fracture appears for a larger displacement than the experimental fracture displacement; Right: numerical fracture appears for a smaller displacement than the experimental fracture displacement.

5. 4 3 2.79 10 

 4310 Parameter Identification using Inverse Method 166 the tolerance of objective function. We set the tolerance of objective function as 6 TolFun 10   . If one of these two criteria is satisfied, the optimization procedure stops. The obtained set of parameters with minimum objective function value are the identified material parameters. In the ductile damage model of Section 5.2, there are Swift hardening law's parameters ( A , 0  and n ), i.e., the strain hardening parameters, and damage parameters ( p D  , r , s and C D ) need to be identified. Since in Lemaitre's damage model 1.0 s  [91], here we set 1.0 s  . Thus, in this paper, we focus on the identification of strain hardening parameters, A , 0  and n , and ductile damage parameters p D  , r and C D . In order to simplify the identification process, these two kinds of parameters (strain hardening parameters and damage parameters) are identified separately. Specifically, the strain hardening parameters ( A , 0  and n ) are identified firstly for the elastoplastic model with Swift hardening law and von Mises yield criterion; then the damage parameters ( p D  , r and C D ) are identified for the ductile damage model. The parameter identification is carried out according to the standard tensile test of A 2017-T4. The aluminum 2017-T4 has the density of kg/m 3 . Its elastic modulus and

5 .

 5 Applications in Material Parameter Identification 167 specimen. The left side is fixed in X direction, and the displacement load is applied at the right side.

Figure 5 . 5

 55 Figure 5.5 Geometry dimension of the test specimen

3 .

 3 The convergence graphs and iterative history of the parameters are illustrated in Figure 5.7. The force-displacement 5.4 Parameter Identification using Inverse Method 168 curves corresponding to identified strain hardening parameters and experimental curve are shown in Figure 5.8.

Figure 5 . 7 Figure 5 . 8

 5758 Figure 5.7 Convergence and iteration graphs of strain hardening parameter identification

Figure 5 . 9 .

 59 Figure 5.9. The force-displacement curves corresponding to identified elastic-plastic damage parameters and experimental curve are shown in Figure 5.10.

Figure 5 . 9 Figure 5 .

 595 Figure 5.9 Convergence and iteration graphs of damage parameter identification

Figure 5 .

 5 Figure 5.11 Comparison of force-displacement curves of elastic-plastic and elastic-plastic damage models with the identified parameters.

6. 2

 2 Organisation de la Thèse 6.2.1 Premier Chapitre: État de l'Art des Techniques d'Optimisation Ce chapitre pré sente un bref aperç u des techniques d'optimisation et de l'é tat de l'art de l'algorithme de l'é volution assisté e par modè le substitutif. On commence par introduire les concepts d'optimisation. Ensuite, on donne un aperç u des algorithmes d'optimisation, y compris les algorithmes avec dé rivé es et ceux sans dé rivé es. Enfin, nous pré sentons l'algorithme de l'é volution assisté e par modè le substitutif. Nous examinons le dé veloppement d'algorithmes é volutifs (EAs) assisté s par substitution. Dans l'optimisation é volutive assisté e par substitution, la fonction objectif est remplacé e par des modè les de substitution. Les mé canismes d'incorporation des substituts dans les EAs peuvent ê tre divisé s en mé thodes de remplacement direct et indirect, c'est-à -dire les styles directs et indirects, comme le montre la Figure 6.1. Cette figure montre le panorama de l'optimisation é volutionnaire assisté e par modè le substitutif.

Figure 6 . 1

 61 Figure 6.1 Les mé canismes d'incorporation des modè les de substitution dans les EAs

15 :

 15 Test d'arrêt vérifié Parmi les ES, la straté gie d'é volution avec l'adaptation de matrice de covariance (CMA-ES) est une straté gie d'é volution trè s dé veloppé e et est devenue une norme pour l'optimisation é volutive continue. Il s'agit d'un puissant algorithme d'optimisation et se comporte particuliè rement bien dans des problè mes de back-box multimodaux non lisses. Le CMA-ES adopte l'opé rateur de mutation corré lé e, ce qui en fait un algorithme de haut niveau comparé à d'autres algorithmes qui utilisent la mutation isotrope. Dans CMA-ES, deux techniques, à savoir l'adaptation de la matrice de covariance (CMA) et l'adaptation cumulative des é chelons (CSA), sont utilisé es pour adapter respectivement la matrice de covariance de la mutation et la taille de l'é chelon. Le CMA-ES est choisi comme la straté gie d'é volution de la base pour les ES assisté es par substitution dans ce travail, en raison de sa puissance et de son succè s dans l'optimisation continue avec des fonctions objectifs implicites. L'algorithme CMA-ES est donné dans l'Algorithme 6.1.

Figure 6 . 2

 62 Figure 6.2 Processus de modé lisation de substitution

6. 2 . 4 . 2 Figure 6 . 3 6. 2 . 4 . 3 Figure 6 . 4 6. 2 . 4 . 4 6. 2 . 4 . 5

 2426324364244245 Figure 6.3 Illustration de la pré sé lection dans la straté gie d'é volution assisté e par modè le de substitution:(1) la norme (,)-ES et (2) l'ES assisté e par la pré -sé lection.

4 )6. 2 . 4 . 7

 4247 KA-CMA-ES utilisant le contrôle basé sur la gé né ration, qui comprend KA-CMA-ES utilisant FGC (contrôle basé sur la gé né ration fixe) et AGC (contrôle adaptatif basé sur la gé né ration).6.2.4.6 ÉtudesExpé rimentales de KA-CMA-ES En utilisant diffé rents modè les d'exploitation dans KA-CMA-ES, diffé rents algorithmes ont é té dé veloppé s. Ceux é tudié s dans cette thè se peuvent ê tre divisé s en quatre groupes : KA-CMA-ES utilisant la pré sé lection qui comprend KA-CMA-ES utilisant PS et CPS, KA-CMA-ES en utilisant un contrôle individuel qui comprend KA-CMA-ES utilisant FIC et MIC, KA-CMA-ES en utilisant une procé dure approximative de classement, et KA-CMA-ES utilisant un contrôle basé sur la gé né ration qui comprend KA-CMA-ES utilisant FGC et AGC. Dans les é tudes expé rimentales de KA-CMA-ES, diffé rents aspects des algorithmes KA-CMA-ES ont é té examiné s en exé cutant les expé riences sur 40 problè mes de benchmark (dont 12 fonctions de test avec des dimensions diffé rentes). Le taux de succè s (SR), la performance de succè s (SP) et/ou la performance d'accé lé ration (SPU) des expé riences sont trois mesures pour é valuer et analyser la performance des algorithmes. Ré sultats Expé rimentaux Dans l'é tude de KA-CMA-ES en utilisant la pré sé lection sans contrôle d'impact de modè le (PS), diffé rentes mé thodes de sé lection d'ensembles d'entraî nement, y compris 'Ré cemment' (les points ré cemment é valué s), 'kNN' Moyenne basé e sur la distance de Mahalanobis), et 'Interval' (la mé thode proposé e d'intervalle de confiance), ont é té é tudié s. Les ré sultats sont montré s dans la Figure 6.5. Les ré sultats de cette é tude ont dé montré que la mé thode d'intervalle de confiance proposé e pour la sé lection d'ensembles d'apprentissage semble supé rieure aux mé thodes « Ré cemment» et « kNN» pour la sé lection d'ensembles d'apprentissage. KA-CMA-ES utilisant PS avec 'Interval' a un taux de ré ussite plus é levé et une meilleure performance de succè s que celle utilisant 'Ré cemment' et 'kNN'. Ainsi, la mé thode de l'intervalle de confiance est suggé ré e pour la sé lection des ensembles d'entraî nement dans la straté gie d'é volution assisté e.

Figure 6 . 5

 65 Figure 6.5 Taux de ré ussite moyen et des performances d'accé lé ration de KA-CMA-ES en utilisant la pré sé lection sans contrôle de l'impact du modè le (PS).

Figure 6 . 6

 66 Figure 6.6 Taux moyen de ré ussite et taux de pré sé lection supé rieur à celui du contrôle de l'impact du modè le à l'aide de diffé rentes mesures de la qualité du modè le.

Figure 6 . 7

 67 Figure 6.7 Taux moyen de ré ussite et accé lé ration de la pré sé lection sans et avec le contrôle de l'impact du modè le.

Figure 6 . 8 Figure 6 . 9

 6869 Figure 6.8 Taux moyen de ré ussite (SR) et performance d'accé lé ration (SPU) de KA-CMA-ES en utilisant un contrôle individuel fixe avec diffé rentes mesures.

Figure 6 .

 6 Figure 6.10 Moyenne SR et SPU de KA-CMA-ES en utilisant un contrôle individuel fixe et mixte (FIC et MIC).

Figure 6 .

 6 Figure 6.12 Taux moyen de ré ussite et accé lé ration de KA-CMA-ES en utilisant FGC et AGC.

6. 2

 2 Figure 6.13 Comparaison du taux de ré ussite moyen (SR) des algorithmes KA-CMA-ES.

Figure 6 . 2 6. 2 . 5 . 3

 62253 Figure 6.14 Comparaison de la performance d'accé lé ration moyenne (SPU) des algorithmes KA-CMA-ES.

Figure 6 .Figure 6 . 16 Figure 6

 66166 Figure 6.15 La procé dure d'identification des paramè tres

Table 4 .

 4

11 Results of KA-CMA-ES using Mixed Individual-based Control (MIC). .................... 126 Table 4.12 Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using mixed individual-based control. ....................................................................................................... 127 Table 4.13 Success rate (SR) and speedup performance (SPU) comparison of KA-CMA-ES using Fixed Individual-based Control (FIC) and Mixed Individual-based Control (MIC). ............ 129

Table 4 .

 4 14 Average success rate and speedup performance of KA-CMA-ES using Fixed and Mixed Individual-based Control (FIC and MIC) .............................................................................. 130

Table 4 .

 4 15 Results of KA-CMA-ES using Approximate Ranking Procedure (ARP) .................... 132

Table 4 .

 4 16 Average success rate and speedup of KA-CMA-ES using ARP. ................................. 133

Table 4 .

 4 17 Results of KA-CMA-ES using FGC and AGC ............................................................ 134

Table 4 .

 4 [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF] Average success rate and speedup of KA-CMA-ES using FGC and AGC. ................ 135

Table 4 .

 4 [START_REF] Jin | Surrogate-assisted evolutionary computation: Recent advances and future challenges[END_REF] Performance comparison of KA-CMA-ES algorithms, including PS-Interval,CPS-

w Q , rank CPS- , FIC-EI, MIC-SLB, ARP-EI, AGC- w Q and rank AGC- . ....................................... 137 XVIII

Table 4 .

 4 20 Comparison of Average success rate (SR) and speedup performance (SPU) of KA-CMA- . .........................................................................................................

	ES algorithms, including PS-Interval, CPS-w Q ,	CPS-	rank	, FIC-EI, MIC-SLB, ARP-EI,
	AGC-

w Q and rank

AGC-

Table 5 .

 5 1 Equations of elastic-plastic damage model ....................................................................

Table 5 .

 5 2 Stress update algorithm for elastic-plastic damage model .............................................

Table 5 .

 5 [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF] Results and computational costs of strain hardening parameter identification ..............

Table 5 .

 5 [START_REF] Koziel | Computational Optimization[END_REF] Results and computational costs of damage parameter identification ............................ Table 6.1 Equations constitutives du modè le d'endommagement élastique-plastique .................. Table 6.2 Ré sultats et coûts de calcul de l'identification des paramè tres de durcissement ............ Table 6.3 Résultats et coûts de calcul de l'identification des paramètres d'endommagement .......Algorithm 1.1 Steepest Descent ....................................................................................................... Algorithm 1.2 Newton's method ...................................................................................................... Algorithm 1.3 Levenberg-Marquardt (LM) method ........................................................................ Algorithm 1.4 Nelder-Mead Simplex Algorithm ............................................................................. Algorithm 1.5 Trust-Region Algorithm ........................................................................................... Algorithm 1.6 General Framework of an Evolutionary Algorithm ..................................................

	Algorithm 2.1 The   1 1 -ES 
	Algorithm 2.2 The 	 , -ES
	Algorithm 2.3 The 	 , -CMA-ES

XIX

List of Algorithms

with 1/5th Success Rule .................................................................... I    with Cumulative Step-Size Adaptation ..................................... w    ..........................................................................................

3.3 Evolutionary Algorithms

	Evolutionary algorithms (EAs) are a class of stochastic derivative-free optimization
	methods. EAs use mechanisms inspired by biological evolution, such as reproduction,
	mutation, recombination, and selection. Among EAs, genetic algorithms (GAs), evolution

strategies (ESs) and differential evolution (DE) are popular and commonly applied in research and engineering problems.

1.2 Overview of Optimization Algorithms

2. Evolution Strategies

	This chapter gives a comprehensive description of evolution strategies (ESs). A short
	introduction of ESs is provided in Section 2.1. Then the main principles and evolutionary
	operators used in ESs are introduced in Section 2.2. Subsequently, Section 2.3 presents
	parameter control of ESs with three ES algorithms, (1+1)-ES, 	 , -ES    I	with
	Cumulative Step-Size Adaptation, and Covariance Matrix Adaptation Evolution Strategy
	(CMA-ES).

  decay factor c  is used to decrease the importance of previously performed steps with time, and the factor

						 cc    2 	is set to normalize the variance of
	 p such that   1 g   p	d	  , 0 I . Finally, the step-size  is updated by using the evolution path
	 p as expressed in Line 11. The step-size increases if the length of the evolution path	p	  1 g  
	is longer than the expected length of the evolution path under random selection
		  , d 0I , and decreases otherwise. The expectation of 	  , d 0I , i.e.		  , d 0I , 
	is approximated by	 1 d  11 4 d 21 d	2		. The damping factor d  is used to control the change of
	the step-size.
	Algorithm 2.2 The 	 , -ES    I	with Cumulative Step-Size Adaptation
	1:	given:

2.3.3 Covariance Matrix Adaptation (CMA)

	algorithms that use isotropic mutation. In CMA-ES, two techniques, namely the covariance
	matrix adaptation (CMA) and the cumulative step-size adaptation (CSA), are used for
	adapting the covariance matrix of mutation and the step-size, respectively.
	The 	 , -CMA-ES    w	is detailed in the remaining part of this subsection. Briefly
	speaking, in generation loop of the 	 , -CMA-ES    w	(simply denoted as CMA-ES
	hereafter),  offspring are generated by mutation based on a single solution vector (the mean
	vector of mutation distribution), comma selection is then performed to select  (   ) best
	individuals according to their fitness, the weighted recombination of selected individuals
	creates the single new solution. The generation loop stops until the termination criterion is
	fulfilled. This generation loop of CMA-ES is presented in details in the following.
				intermediate recombination
	10:		  11         : 1 12 g g g i i c c c                 p p z m	//update evolution path
	11:		  1 gg            1 exp , g c d       1 p       0I 	//step-size update
	12:		gg 	1
	13:	until termination criterion is fulfilled
		The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [39, 40], is a highly
	developed evolution strategy and has become a standard for continuous black-box
	evolutionary optimization. It is a powerful optimization algorithm and performs especially
	well in non-smooth, multimodal back-box problems. The CMA-ES adopts the correlated
	mutation operator, which makes it a high-level algorithm compared with previous described

  comprehensive study can be found in[58, 59].

	Then, training a network is to estimate the unknown parameters in (3.32). Usually,
	MLP utilizes a supervised learning technique called backpropagation for training the
	network [2, 59]. One of the drawbacks associated with neural networks for function

approximation is the fairly large number of parameters that need to be prescribed by the user, thereby demanding adequate user experience in implementing MLP. These prescribed 3. Surrogate Modeling of Computer Experiments 63

   is the bias. A small w in Equation(3.34) means that the regression is flat[START_REF] Smola | A tutorial on support vector regression[END_REF]. One way to ensure the flatness is to minimize the norm,

	  fb    T  x w x	(3.34)
	11 , yy , , , nn x x is a map (or transform) that transforms the input space (when ,  X y x , y  a general SVR model can be given as where    to some feature space , i.e., map   :  x ,  w is the column vector of weights, where d  x and d  x 2 T ,  w w w w w where ,  denotes the dot product. Considering minimizing this norm under the restriction of deviation, the , is d ) construction of the SVR model is reduced to the following convex optimization problem:

b

  3.4 Model Validation/Model PerformanceAssessment 70 sufficient to measure surrogate model accuracy, and the validation with an additional data set is hence recommended[START_REF] Lin | An Efficient Robust Concept Exploration Method and Sequential Exploratory Experimental Design[END_REF].CV is an extremely popular method for verifying the prediction capability of a surrogate model, when additional validation points are unavailable or expensive. The disadvantage of this method is that the surrogate model has to be constructed more than once. For example, in a leave-one-out CV with n samples, the surrogate model need to be trained for n times. If the construction of the surrogate model is expensive, CV would be impractical. This challenge has stimulated the advent and development of surrogate-assisted ES ( also called fitness approximation in ES), which becomes a promising solution to reduce the computation cost of ES. In surrogate-assisted evolution strategies, generally, with the evaluated points, a surrogate model is trained to approximate the real (or original) fitness function and then used together with the real fitness function to guide the search of promising solutions. In this way, the computational cost is reduced because the evaluation of a surrogate is much cheaper than that of the expensive fitness function.

	Among the ESs, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is a
	highly developed evolution strategy and has become a standard for continuous black-box
	evolution strategy. Consequently, this chapter concentrates Kriging-Assisted Covariance
	Matrix Adaptation Evolution Strategy (CMA-ES), which is abbreviated as KA-CMA-ES in
	the thesis.

evolutionary optimization. It is a powerful optimization algorithm and performs especially well in non-smooth, multimodal back-box problems. The CMA-ES adopts the correlated mutation operator, which makes it a high-level algorithm compared with other algorithms that use isotropic mutation. In CMA-ES, two techniques, namely the covariance matrix adaptation (CMA) and the cumulative step-size adaptation (CSA), are used for adapting the covariance matrix of mutation and the step-size, respectively. The CMA-ES is selected as the basis evolution strategy for surrogate-assisted ES in this work owing to its powerfulness and success in continuous black-box optimization.

Compared with other surrogate models like RBF, MLP model, the advantage of Kriging model is that it not only predicts fitness value but also provides the standard deviation for the predicted fitness without additional computational cost. The standard deviation of predicted fitness indicates the uncertainty in fitness approximation. This information is valuable for evolution control and improving the quality of surrogate in evolutionary computation process. That is the reason for selecting Kriging model for surrogate-assisted

  essential goal of surrogate-assisted evolutionary optimization is to use one or several surrogate models f of the fitness function (objective function) f to improve the quality of the search especially in terms of number of functions evaluations required to reach the optimum. It is usually assumed that the evaluation of f is expensive in terms of time or money and the learning (or training) of f is relatively cheap. Surrogate model f is trained

	(learned/built) from a set of pairs		  xx which is known as training data set or simply  , f
	training set. Thus, the selection of training set has influence on the quality of the surrogate

model developed thereof. Generally speaking, a surrogate model with good quality ensures good performance in surrogate-assisted optimization. In this section, a short overview of existed works on surrogate-assisted evolution strategies (ES) is presented and some challenges in this field are described subsequently.

  -selection (pre-screening) of promising solutions (based on approximate model) is another popular strategy of exploiting information from the surrogate model f in Evolution Strategies. In pre-selection strategy, Pre   (  is the population size of ES) individuals are generated each generation, then all Pre  individuals are evaluated by approximate model f , after that,  out of Pre  best individuals are selected to evaluate by the original fitness function.

  1, evaluates the quality of the surrogate model by its consistency in ranking the population rather than its statistical accuracy. The approximate ranking procedure in Algorithm 4.1 can be seen as an adaptive evolution control mechanism to determine the number of controlled individuals in every generation as follows: individuals are successively selected to be evaluated according to their approximate fitness and then are added to the training set until the surrogate-based selection of the parents remains unchanged in two iteration cycles. Then, based on approximate ranking procedure and using locally weighted regression as surrogate, the local meta-model CMA-ES (lmm-CMA-ES) was firstly proposed by Kern et al.[START_REF] Kern | Local Meta-models for Optimization Using Evolution Strategies[END_REF] and later extended for large population size by Bouzarkouna et al.[START_REF] Bouzarkouna | Investigating the Local-Meta-Model CMA-ES for Large Population Sizes[END_REF]. The lmm-CMA-ES is a carefully designed surrogate-assisted CMA-ES algorithm.

	Algorithm 4.1 Approximate Ranking Procedure
	1: given:   k k  xx   1 , ft , , 	(archive which save all evaluated data points).
	2: approximate: build surrogate models f based on , and predict	f	kk   ˆ, f  x k	1, ,	
	3: rank and determine the parent set

  Even though surrogate-assisted ES has achieved considerable improvements over the past decade, there are still many open issues need to be addressed. These open issues mainly lies in i) surrogate model learning, i.e. how to use the previously evaluated data and training the surrogate model, ii) model quality measure, i.e. how to evaluate the quality of the surrogate model, and iii) model exploitation, i.e., how to exploit the built surrogate model. There are several open questions for strategies of surrogate model exploitation, i.e., pre-selection and evolution control. The following will briefly describe some open issues of training set selection and model exploitation, which will be investigated and addressed in this work in One can use all the previously evaluated points as training set to construct surrogate models. Generally, for the same problem and using same modeling techniques, the size of training set is larger, the performance of the trained model is better. A model learned from all previously points may have better performance than one learned from just a part of existed points. However, when all the evaluated points are taken as training data, the size of training set increases during the evolutionary search, and this require higher computational cost for model training. In surrogate-assisted ES, usually, surrogate model need to retrained for large number of times, so the high training cost of the surrogate model is not desirable.

	subsequent sections.
	4.2.2.1 Training Set Selection
	Let training

Surrogate model learning is an essential step for surrogate-assisted ES. In surrogate model learning, the key task is learning one or several appropriate surrogate models from the evaluated candidates. Training data set selection has influence on the quality and the training cost of surrogate model f (besides, the modeling method also has significant on the training cost of f ). Thus, it is an important issue in surrogate model learning.

Another issue that relate to both design and use of surrogate models for evolutionary computation is their quality. A surrogate with poor quality may guides the search to false optimum and thus, model quality assessment is of importance.

After surrogate model has been trained, the key aspect of surrogate-assisted ES is to exploit information from the surrogate model f . The two main strategies of surrogate model exploitation applied in surrogate-assisted ES are pre-selection and evolution control. n denotes the number of data points in training set training . Three popular ways to select the training set (points/individuals for model training) are: 1. All Previously Evaluated Points 2. Recently Evaluated Points A simple but efficient strategy is to choose training n the most recently evaluated points as training data set. This is a reasonable strategy because that the recently evaluated points are relatively close to the actual search subspace. However, in some cases, this strategy may be insufficiently "greedy" because it may "forgets" promising individuals generated more than training n

  [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF] is used as a normalization in Equation (4.2). If all  parents individuals are selected correctly,

	the measure select  reaches its maximum value of select 1   . The negative values of select 

indicate that the selection based on the surrogate (model-based selection) is worse than a random selection. The measure select  only evaluates the absolute number of correctly selected individuals. However, this measure does not consider the rank of the selected individuals. In case of select 1   , the measure select  does not indicate whether the ()-th best individual or the worst individual has been selected, which may have significant influence on the evolution process. Therefore, the measure select  is extended to include the rank of the selected individuals that are determined based on the original fitness function values. In the extended measure select  , the surrogate model gets a grade of m   , if the m-th best individual based on the original fitness function is selected. Thus, the quality of the surrogate model can be indicated by summing up the grades of the selected individuals, which is demoted by . Apparently, if all  individuals are selected correctly, reaches its maximum:

Nearest Neighbor Points to Distribution Mean based on Mahalanobis Distance

  Training set selection affects both the quality of surrogate model learned from it and the cost of model training. Thus, it is a subject worth studying. In previous section, selecting training set as all previously evaluated points, the recently evaluated points, and nearest neighbor points are already described. In this work, we propose a new approach for training set selection based on properties of multivariate normal random variables. Keep in mind that our goal of learning surrogate models is to predict fitness values of current offspring population. Therefore, it is rational to choose points that are relevant to current offspring population or current search subspace. The previously mentioned Recently Evaluated Points and k-Nearest Neighbor Points for training set selection are based on this 4.3 Address Some Open Issues of Surrogate-Assisted ES 84principle. However, in Recently Evaluated Points approach, no measurement is used to assess the relevance between recently evaluated points and current search subspace. So, it cannot guarantee that the selected recently evaluated points all are relevant to current search subspace. In selection of training set by the k-Nearest Neighbor Points, the Euclidean distance is usually used as a measurement to indicate the similarity between evaluated points and the query point in such a way that previously evaluated points with smaller Euclidean distance to the query point are chosen. For evolution strategy algorithms using isotropic mutation (uncorrelated mutation), using Euclidean distance as similarity measure metric is appropriate. When correlated mutation operator is used in evolution strategy algorithms, such as the CMA-ES, the Euclidean distance may not be useful on correlated data since there is no adjustment for the covariance. In this situation, the Mahalanobis distance becomes a more proper similarity measure metric. Therefore, we firstly suggest to use the Mahalanobis

distance instead of Euclidean distance in k-Nearest Neighbor Points approach for training set selection and another approach of using confidence interval is proposed. 4.3.1.1 k-Without loss of generality, we consider the   , -CMA-ES w    (this algorithm has presented in Section 2.3.3) which uses a general (or correlated) mutation operator. In the evolution loop of CMA-ES, for generation g ,  offspring are generated by mutation as:

Model Impact Control based on Different Model Quality Measures

					4. Kriging-Assisted CMA Evolution Strategy
					  1 Pre g  	       	    Pre Pre gg ww       rand Pre max rand rand Pre rand , if , otherwise g ww rand ww g g ww w QQ QQ QQ QQ Q       	(4.15)
	where Pre   is the adaptation rate, which is set by the user.
	In pre-selection with model impact control, different model quality measures can be
	used. In our work, three model quality measures are included. These three model quality
	measures are: 1) proposed selection-based model quality measure using recombination
	weights	w Q ; 2) selection-based model quality measure selection Q	proposed by Ulmer et al. [84];
	3) rank correlation coefficient rank
					g QQ 	TS	; and Pre  is decreased
	if   g QQ 	TS	. The uniform formulation of model impact control based on model quality
	measure Q , which includes	w Q , selection Q	and rank  , can be given as:
					  Pre gg   g TS   Pre max TS , if QQ QQ 	TS
	  g w Q is computed and compared with the expected quality of the QQ Q (given in (4.14)). random selection rand w   Pre Pre TS min , otherwise g QQ      1 Pre g TS g QQ  
	If   g ww rand QQ 	the model based selection is better than a random selection, and Pre  should
	be increased. On the contrary, for   ww rand g QQ 	the Pre  should be decreased. When   ww rand g QQ 
	the value of Pre  is kept. Therefore, Pre  is adapted similarly as (4.12), i.e.,
					91

 proposed by Jin et al.

[START_REF] Jin | Quality measures for approximate models in evolutionary computation[END_REF]

. The strategy of updating the pre-selection population size Pre  is expressed by Equation (4.12) or (4.15). In the model impact control rule, three critical values of model quality measures and the model quality evaluated at current generation are used. These three critical values are the minimum, maximum and a threshold value of model quality, which are respectively noted by min Q , max Q and TS Q . The size of pre-selection population Pre  is increased when the model quality of current generation is larger than the threshold value, i.e.,  

Table 4 .

 4 1 Critical values of model quality measures Model quality measure Q Minimum value min Q

  In this work, we comprehensively investigate the criteria which based on the model prediction   f x or/and estimated standard deviation  

	ŝ x to select
	individuals for re-evaluation. The criteria used to select individuals which will be controlled
	are called metrics hereafter in this thesis. Five different metrics, which have been used in
	surrogate-based optimization [54, 57] but relatively new in surrogate-assisted evolution
	computation, including the mean of model prediction (Mean), estimated standard deviation
	(SD) of prediction, statistical lower bound (SLB), probability of improvement (POI), and
	expected improvement (EI), are introduced and studied in KA-CMA-ES using fixed
	individual-based control. These five metrics are defined by
	individual-based control (FIC), usually, the individuals with best predicted fitness
	are selected to be controlled. Compared with other surrogate models like RBF, MLP model,

the advantage of Kriging model (or Gaussian process model) is that it not only predicts fitness value but also provides the standard deviation for the predicted fitness without additional computational cost.

  , MIC) which combines both features of fixed and adaptive control in some way, is proposed and presented below. Similarly to fixed individual-based control, in this mixed individualbased evolution control strategy,   individuals are chosen by a best strategy for re- evaluation in each generation. However, the value of  is vary between the user defined lower

	4.3 Address Some Open Issues of Surrogate-Assisted ES
	fail of converging to the global optimum and too large min  would require more exact fitness
	function evaluations.
	bound min  and upper bound max  , i.e.	min 01  max       , which is mainly different from
	the FIC. In each generation, min   individuals selected by a best strategy are unconditionally
	re-evaluated using original fitness function  . f x Up to  	max			min			more offspring which
	take on better fitness values (estimated by surrogates   f x ) than the current best solution
	can be re-evaluated too. With these newly evaluated points, the surrogate model is retrained
	(updated) and the fitness values of uncontrolled individuals are re-estimated by the updated
	kk x surrogate model. k   z	m	ggg C  t f x	C	x
	2: training set selection:	T		training_set_selection	 
	3: model training:	ˆmodel_training f 	  T
	4:	for	k			1, ,		do
	5:	C	kk   C  x	computing by model	metric f //compute the metric values of individuals
	6:	end for
	7:	controlled	   1: 2: uncontrolled : , , , ,         x x x	 x	    	1:  	, ,	x	:  		//best strategy based on   C x
	8:	for	x	k 	controlled	do
	9:		  kk ff  x		//evaluated by original fitness function
	10:	1 Above described mixed individual-based control strategy can be regarded as a tt 
	11: combination of fixed control with control frequency of min   , kk f  x   and an adaptive control that
	12: controls up to   TT   x max , kk f 	 	min			offspring individuals that has predicted fitness better than
	13: current best solution. In addition, after the re-evaluation of controlled individuals, on-line end for
	14: model updating is performed by retraining the surrogate model based on training set that   T ˆtrain_model f  //update surrogate model
	15:	for	x	k 	uncontrolled	do
	16: adds these newly evaluated points. This could improve the quality of surrogate model and is   kk ff  x //evaluated by surrogate model
	end for good to the evolutionary computation. 17:
	18:	output: Above proposed mixed individual-based control (MIC) is implemented here and given   1 , ,  , , k k k k tf 
	in Algorithm 4.4. The metric Mean and SLB can be used in MIC. The lower bound min  and
	upper bound max  are set as min 0.4   and max 0.8   . This means that at least 40% and up to
	80% individuals are controlled in each generation. Too small min  may lead to the algorithm
												98

xz 4.3.3.2 Mixed Individual-based Control (MIC) With the efforts of drawing the advantages of both fixed and adaptive individual-based evolution control, an individual-based control strategy called mixed Individual-based control (Furthermore, since Kriging model is used as the surrogate model, a metric function   C x that is based on both the model prediction   f x and estimated standard deviation   ŝ x to identify the most promising individuals whose fitness are probably better than the current best fitness. Evidently, besides the mean of prediction (Mean), the statistical lower bound (SLB) can be used in above described mixed individual-based control without any other change of the control method.

Algorithm 4.4 Mixed Individual-based Control (MIC) 1: given:

  best strategy using  

										C x
	8: for	k	   min     	 1 , ,	max    	  do
	9:	if : k Cf  	min	then	// min f is the minimum fitness in
	10:				controlled	 kk  : uncontrolled ,  x controlled x		1: 	, ,	x	:  	
	11:	end if
	12: end for
	13: for	x	k 	controlled	do
	14:		  kk ff  x	//evaluated by original fitness function
	15:	tt 	1
	16:							x	, kk f	
	17:			TT 			x	, kk f	
	18: end for
	19:	ˆtrain_model f 	  T	//update surrogate model
	20: for	x	k 		uncontrolled	do
	21:		  kk ff  x	//evaluated by surrogate model
	22: end for
	23: output:	 tf , , , , k k xz	k	 1 k  
	4.3.4 Modified Approximate Ranking Procedure (ARP)

  . Because it is very difficult to estimate the global fidelity of the surrogate model, the authors proposed to use a local estimation of the model fidelity. Since the evolution strategy generally proceeds with small steps, the current error can be used to estimate the local fidelity of the model and then to Ek is the current model error estimation, k denotes the k-th evolution control cycle of p generations. The current model error, which is estimated before the next control cycle, is measured as follows:

	generations need to be controlled. The  , which indicates the control frequency, is adapted
	by
		  1,     min max min max Ek k E          	(4.22)
	where x   denotes the largest integer that is smaller than x , max  is the maximal value of  ,
	max   , and the minimum min p  usually equals 1 so that the information on the model fidelity
	is always available, max E is the allowed maximal model error and  

determine the frequency at which the original fitness function is used and the surrogate model is updated. Consider there are p generations within an evolution control cycle, and 

4 Kriging-Assisted CMA-ES (KA-CMA-ES) Algorithms

  next generation is control or not. The quality of the surrogate model Q is estimated

	In this section, these issues described in previous section including training set selection,
	pre-selection with model impact control (CPS), fixed individual-based control using metrics
	(FIC), mixed individual-based control (MIC), modified approximate ranking procedure
	(ARP) and adaptive generation-based control (AGC) are incorporated into CMA-ES to
	compose several different KA-CMA-ES algorithms. In training and using the surrogate
	model, variable transformation of correlated variables is performed firstly. Before the
	evolution loop of KA-CMA-ES, a number of solutions (candidates or points) are initially
	in every controlled generation, if	QQ 	TS	( TS Q is the critical value or threshold of model
	quality), the next generation is evaluated by surrogate model; otherwise, the next generation
	is controlled and evaluated by original fitness function. The model quality measurements
	can be the proposed selection-based model quality measure using recombination weights	w Q ,
	selection-based model quality measure selection Q	proposed by Ulmer et al. [84], and the rank
	correlation coefficient rank

 . The threshold of model quality measurements have been given in Table

4

.1. 4.4 Kriging-Assisted CMA-ES (KA-CMA-ES) Algorithms 102 4.

sampled such that surrogate model can be built at the beginning of evolution loop. At the same time, the best candidate in initial sampling are chosen as the start point for evolution loop in KA-CMA-ES rather than set the start point randomly. All these are described in the remaining part of this section.

Table 4

 4

	.2,

Table 4

 4

				.2 Test functions for experimental studies
	Name	Function	Search Space
	Sphere	1	 	2
			1

Pre-Selection 4.5.2.1 Pre-Selection without Model Impact Control (PS)

  which method is used in subsequent experiments of other KA-CMA-ES algorithms. Additionally, the performance of PS algorithms will be compared with the preselection with model impact control (CPS) which are studied in next subsection. Table4.3. In this table, PS-Recently indicates the pre-selection using Recently Evaluated Points method for training set selection, PS-kNN stands for the pre-selection where training set is selected as k-Nearest Neighbor Points to distribution mean based on Mahalanobis distance, and PS-Interval signifies pre-selection using the proposed Confidence Interval method as training set selection.

	4.5 Experimental Studies
	thus to determine The experimental results of KA-CMA-ES using pre-selection without model impact
	control (including success rate (SR), success performance (SP) and speedup performance
	(SPU)), in which training set selection methods 'Recently', 'kNN' and 'Interval' are adopted,
	are listed in
	Neighbor Points to distribution mean based on Mahalanobis distance and Confidence Interval, have been investigated in KA-CMA-ES using pre-selection without model impact control (PS). The size of training set is set as T 2   for Recently Evaluated Points and k-Nearest Neighbor Points methods, and training 99.73%  for Confidence Interval methods. The size of pre-selection population is set as Pre 2   and keeps as constant during the computation. For simplification reason, we use 'Recently' to represents the Recently Three training set selection methods, including Recently Evaluated Points, k-Nearest n p Points to distribution mean based on Mahalanobis distance, and 'Interval' to signify the proposed Confidence Interval method. The goal of running KA-CMA-ES using PS (KA-Evaluated Points method of training set selection, 'kNN' to stand for the k-Nearest Neighbor 114

CMA-ES using pre-selection without model impact control) with different training set selection methods is to investigate the performance of training set selection methods and From Table 4.3, it can be found that, for large majority of the test problems, the speedup performance (SPU) of PS-Interval is larger than that of PS-Recently and PS-kNN, which indicates better performance of PS-Interval. Additionally, PS-Interval can get higher success rate (SR) values on difficult problems such as 11 f and 12 f . It can be initially concluded that the proposed Confidence Interval method for training set selection, generally, work better than 'Recently' and 'kNN' methods.

Table 4 .

 4 [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF] Results of KA-CMA-ES using pre-selection without model impact control (PS)

	Function D	λ	PS-Recently SR SP SPU SR	PS-kNN SP SPU SR	PS-Interval SP SPU SR CMA-ES SP
		2	6	253 1.31 1	255 1.30 1	257 1.29 1	331
	f1	5 10 10 8	718 1.24 1 1535 1.16 1 1521 1.17 1 708 1.25 1	705 1.26 1 1410 1.27 1	888 1784
		20 12	3110 1.06 1 3090 1.07 1	2902 1.13 1	3292
		2	6	326 1.19 1	323 1.20 1	319 1.21 1	387
	f2	5 10 10 8	971 1.16 1 2632 1.12 1 2422 1.21 1 998 1.13 1	1004 1.12 1 2509 1.17 1	1128 2941
		20 12 0.96 6422 1.16 1 6742 1.11 0.96 6634 1.13 1	7470
		2	6	279 1.28 1	283 1.26 1	274 1.30 1	356
	f3	5 10 10 8	796 1.23 1 1850 1.12 1 1808 1.14 1 801 1.22 1	760 1.29 1 1661 1.25 1	980 2070
		20 12	4416 0.98 1 4223 1.02 1	4194 1.03 1	4327
		2	6	325 1.28 1	336 1.24 1	326 1.28 1	417
	f4	5 10 10 8	979 1.25 1 2535 1.13 1 2530 1.13 1 995 1.23 1	941 1.30 1 2307 1.24 1	1224 2858
		20 12	7172 1.04 1 7114 1.04 1	6825 1.09 1	7431
		2	6	202 1.37 1	219 1.26 1	216 1.28 1	277
	f5	5 10 10 8	878 1.31 1 2866 1.24 0.92 3173 1.12 1 888 1.29 1	912 1.26 1 2794 1.27 1	1149 3556
		20 12 0.72 12036 1.16 0.88 9843 1.42 0.84 10444 1.34 0.76 14014
		2	6	578 1.35 1	574 1.36 1	578 1.35 1	778
	f6	5	8	1745 1.29 1 1699 1.33 1	1599 1.41 1	2253
		10 10	4133 1.23 1 4109 1.24 1	3709 1.37 1	5078
		2	6	470 1.41 1	456 1.45 1	449 1.48 1	663
	f7	5	8	1744 1.42 1 1809 1.36 0.96 1675 1.47 0.96 2468
		10 10 0.88 6503 1.27 0.96 5895 1.40 0.96 5073 1.63 0.8 8248
		2	6	553 1.30 1	560 1.28 1	557 1.29 1	717
	f8	5	8	1464 1.39 0.96 1513 1.35 1	1414 1.44 0.92 2040
		10 10 0.96 3163 1.19 1 3052 1.23 1	2810 1.34 0.96 3753
		2 12	390 1.25 1	381 1.28 1	376 1.29 1	486
	f9	5 16	1141 1.30 1 1137 1.31 1	1137 1.31 0.96 1489
		10 20	2492 1.22 0.96 2596 1.17 1	2440 1.25 1	3049
		2 12	846 1.25 1	852 1.24 0.96 901 1.17 1	1054
	f10	5 16	2354 1.33 1 2399 1.31 1	2362 1.33 0.96 3137
		10 20	5144 1.34 1 5136 1.34 1	4924 1.40 0.92 6893
		2 12	454 1.31 1	461 1.30 1	460 1.30 1	597
	f11	5 16 0.8 1618 1.05 0.92 1405 1.21 1	1299 1.31 0.96 1700
		10 20 0.88 3152 1.41 0.76 3673 1.21 0.92 2949 1.51 0.76 4457
	f12	2 50 0.92 1508 1.57 0.96 1452 1.63 0.92 1496 1.59 0.68 2372 5 140 0.6 9722 1.03 0.76 8338 1.20 0.84 7452 1.34 0.8 10020

Table 4 .

 4 [START_REF] Koziel | Computational Optimization[END_REF] Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using preselection without model impact control (PS).

	Category	Average Success Rate (SR) PS-Recently PS-kNN PS-Interval CMA-ES PS-Recently PS-kNN PS-Interval Average Speedup Performance (SPU)
	D=2	0.993	0.997	0.990	0.973	1.322	1.316	1.318
	D=5	0.950	0.970	0.983	0.963	1.251	1.267	1.321
	D=10	0.975	0.964	0.989	0.949	1.221	1.216	1.335
	D=20	0.936	0.976	0.960	0.952	1.080	1.133	1.145
	Unimodal	0.983	0.991	0.989	0.982	1.221	1.230	1.277
	Multimodal	0.940	0.951	0.974	0.923	1.282	1.290	1.347
	Overall	0.968	0.977	0.984	0.961	1.242	1.251	1.302

Table 4 .

 4 [START_REF] Chong | An introduction to optimization[END_REF] Results of KA-CMA-ES using pre-selection with model impact control (CPS).

	Function D	λ	SR	CPS-w Q SP	SPU SR	selection CPS-Q SP SPU SR	rank SP SPU SR CPS- CMA-ES SP
		2	6	1	206 1.61 1	209 1.58 1	202 1.64 1	331
	f1	5 10 10 8	1 1	541 1.64 1 1122 1.59 1	544 1.63 1 1159 1.54 1 1220 1.46 1 545 1.63 1	888 1784
		20 12	1	2767 1.19 1	2756 1.19 1 2871 1.15 1	3292
		2	6	1	256 1.51 1	264 1.47 1	258 1.50 1	387
	f2	5 10 10 8	1 1	797 1.42 1 1864 1.58 1	753 1.50 1 1829 1.61 1 1782 1.65 1 814 1.39 1	1128 2941
		20 12	1	3778 1.98 1	3974 1.88 1 4330 1.73 1	7470
		2	6	1	220 1.62 1	216 1.65 1	222 1.60 1	356
	f3	5 10 10 8	1 1	589 1.66 1 1342 1.54 1	598 1.64 1 1389 1.49 1 1468 1.41 1 595 1.65 1	980 2070
		20 12	1	5307 0.82 1	5422 0.80 1 4140 1.05 1	4327
		2	6	1	254 1.64 1	256 1.63 1	252 1.65 1	417
	f4	5 10 10 8	1 1	710 1.72 1 2102 1.36 1	719 1.70 1 2084 1.37 1 2224 1.29 1 723 1.69 1	1224 2858
		20 12 0.6 15281 0.49 0.6 14963 0.50 1 6754 1.10 1	7431
		2	6	1	163 1.70 1	165 1.68 1	163 1.70 1	277
	f5	5 10 10 8	1 1	572 2.01 1 1588 2.24 1	594 1.93 1 1693 2.10 1 1712 2.08 1 596 1.93 1	1149 3556
		20 12	1	4279 3.28 1	4548 3.08 1 4852 2.89 0.76 14014
		2	6	1	448 1.74 1	461 1.69 1	466 1.67 1	778
	f6	5	8	1	1236 1.82 1	1291 1.75 1 1358 1.66 1	2253
		10 10	1	2881 1.76 1	2853 1.78 1 3305 1.54 1	5078
		2	6	1	308 2.15 1	311 2.13 1	320 2.07 1	663
	f7	5	8 0.96 1104 2.24 1	1105 2.23 0.96 1222 2.02 0.96 2468
		10 10 0.96 4013 2.06 0.92 4444 1.86 0.96 4555 1.81 0.8 8248
		2	6	1	431 1.66 1	438 1.64 0.96 463 1.55 1	717
	f8	5	8	1	1043 1.96 1	1067 1.91 1 1070 1.91 0.92 2040
		10 10	1	2096 1.79 1	2153 1.74 1 2364 1.59 0.96 3753
		2 12	1	342 1.42 1	341 1.43 1	344 1.41 1	486
	f9	5 16	1	977 1.52 1	985 1.51 1	989 1.51 0.96 1489
		10 20	1	2046 1.49 0.96 2180 1.40 1 2198 1.39 1	3049
		2 12 0.52 1493 0.71 0.4 1878 0.56 0.92 865 1.22 1	1054
	f10	5 16	1	1899 1.65 1	1946 1.61 1 2141 1.47 0.96 3137
		10 20	1	3946 1.75 1	4058 1.70 1 4717 1.46 0.92 6893
		2 12	1	399 1.50 1	404 1.48 1	416 1.44 1	597
	f11	5 16 0.96 1134 1.50 1	1089 1.56 0.92 1204 1.41 0.96 1700
		10 20 0.92 2460 1.81 0.88 2635 1.69 0.8 2970 1.50 0.76 4457
	f12	2 50 0.8 1700 1.40 0.88 1532 1.55 0.96 1417 1.67 0.68 2372 5 140 0.64 9092 1.10 0.8 7376 1.36 0.84 7222 1.39 0.8 10020

Table 4 .

 4 6 Average success rate (SR) and speedup (SPU) performance of KA-CMA-ES using preselection with model impact control (CPS).

	Category	CPS-w Q	Average Success Rate (SR) selection CPS-Q rank CPS-	Average Speedup (SPU) w Q CMA-ES CPS-selection CPS-Q CPS-	rank
	D=2	0.943	0.940	0.987	0.973	1.554	1.540	1.594
	D=5	0.963	0.983	0.977	0.963	1.687	1.695	1.637
	D=10	0.989	0.978	0.978	0.949	1.724	1.662	1.561
	D=20	0.920	0.920	1.000	0.952	1.549	1.490	1.581
	Unimodal	0.982	0.982	0.997	0.982	1.706	1.669	1.651
	Multimodal 0.917	0.923	0.957	0.923	1.518	1.510	1.493
	Overall	0.959	0.961	0.983	0.961	1.640	1.614	1.596

Table 4 .

 4 7 Comparison of KA-CMA-ES using pre-selection without and with model impact control, including PS-Interval, CPS-w Q ,

	4. Kriging-Assisted CMA Evolution Strategy

Table 4 .

 4 [START_REF] Papadrakakis | Survey of quasi-Newton methods with reduced storage[END_REF] Average success rate and speedup performance of KA-CMA-ES using pre-selection without and with model impact control, including PS-Interval, CPS-w Q ,

	selection CPS-Q	and	rank CPS- .

Table 4 .

 4 [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] Results of KA-CMA-ES using Fixed Individual-based Control (FIC) with different metrics.

	4.5 Experimental Studies

Table 4 .

 4 10 Average success rate and speedup performance of KA-CMA-ES using fixed individualbased control (FIC) with different metrics.

	Category	Average Success Rate (SR) FIC-Mean FIC-SD FIC-SLB FIC-POI FIC-EI CMA-ES FIC-Mean FIC-SD FIC-SLB FIC-POI FIC-EI Average Speed (SPU)
	D=2	0.977	0.997	0.990	0.980 0.993 0.973	1.857	1.862	1.902	1.858 1.920
	D=5	0.923	0.937	0.940	0.960 0.953 0.963	1.757	1.587	1.806	1.866 1.824
	D=10	0.815	0.935	0.909	0.876 0.971 0.949	1.456	1.327	1.606	1.554 1.738
	D=20	0.944	0.960	0.920	0.952 0.968 0.952	1.473	1.319	1.455	1.536 1.576
	Unimodal	0.948	0.982	0.972	0.971 0.988 0.982	1.661	1.554	1.701	1.716 1.747
	Multimodal	0.846	0.911	0.891	0.889 0.943 0.923	1.683	1.584	1.800	1.775 1.892
	Overall	0.912	0.957	0.944	0.942 0.972 0.961	1.669	1.564	1.736	1.737 1.798

Table 4 .

 4 [START_REF] Rios | Derivative-free optimization: A review of algorithms and comparison of software implementations[END_REF] Results of KA-CMA-ES using Mixed Individual-based Control (MIC).

	Function	D	λ	SR	MIC-Mean SP	SPU	SR	MIC-SLB SP	SPU	CMA-ES SR SP
		2	6	1	131	2.53	1	131	2.53	1	331
	f1	5 10	8 10	1 1	395 1007	2.25 1.77	1 1	449 1140	1.98 1.56	1 1	888
		20	12	1	2387	1.38	1	2590	1.27	1
		2	6	1	148	2.61	1	147	2.63	1	387
	f2	5 10	8 10	1 1	473 1208	2.38 2.43	1 1	459 1256	2.46 2.34	1 1
		20	12	1	3251	2.30	1	3664	2.04	1
		2	6	1	143	2.49	1	143	2.49	1	356
	f3	5 10	8 10	1 1	428 1318	2.29 1.57	1 1	505 1342	1.94 1.54	1 1	980
		20	12	1	4146	1.04	1	3624	1.19	1
		2	6	1	165	2.53	1	163	2.56	1	417
	f4	5 10	8 10	1 1	535 1882	2.29 1.52	1 1	617 1987	1.98 1.44	1 1
		20	12	1	6982	1.06	1	6258	1.19	1
		2	6	1	114	2.43	1	115	2.41	1	277
	f5	5 10	8 10	1 0.84	500 1889	2.30 1.88	1 0.84	506 1986	2.27 1.79	1 1
		20	12	0.6	8093	1.73	0.6	8140	1.72	0.76	14014
		2	6	0.96	323	2.41	1	340	2.29	1	778
	f6	5	8	1	1812	1.24	1	1270	1.77	1
		10	10	0		0.0	0.96	4176	1.22	1
		2	6	1	265	2.50	1	280	2.37	1	663
	f7	5	8	0.92	1138	2.17	0.92	1251	1.97	0.96
		10	10	0.84	6086	1.36	0.96	4293	1.92	0.8
		2	6	0.88	313	2.29	0.96	308	2.33	1	717
	f8	5	8	0.92	835	2.44	1	842	2.42	0.92
		10	10	0.56	3573	1.05	0.96	2189	1.71	0.96
		2	12	1	218	2.23	1	234	2.08	1	486
	f9	5	16	1	587	2.54	0.96	667	2.23	0.96
		10	20	0.88	1590	1.92	0.92	1654	1.84	1
		2	12	0.64	763	1.38	1	501	2.10	1
	f10	5	16	0.84	1464	2.14	0.96	1349	2.33	0.96
		10	20	0.2	15116	0.46	0.52	5870	1.17	0.92
		2	12	0.96	270	2.21	1	263	2.27	1	597
	f11	5	16	0.84	770	2.21	0.8	840	2.02	0.96
		10	20	0.52	3029	1.47	0.64	2584	1.72	0.76
	f12	2 5	50 140	0.84 0.24	1181 12172	2.01 0.82	0.84 0.4	1165 8568	2.04 1.17	0.68 0.8	10020

Table 4 .

 4 [START_REF] Conn A. R | Recent progress in unconstrained nonlinear optimization without derivatives[END_REF] Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES using mixed individual-based control.

	Category	Average Success Rate (SR) MIC-Mean MIC-SLB CMA-ES	Average Speedup (SPU) MIC-Mean MIC-SLB
	D=2	0.940	0.983	0.973	2.302	2.341
	D=5	0.897	0.920	0.963	2.090	2.046
	D=10	0.713	0.891	0.949	1.449	1.661
	D=20	0.920	0.920	0.952	1.503	1.483
	Unimodal	0.929	0.972	0.982	1.961	1.957
	Multimodal	0.737	0.854	0.923	1.798	1.960
	Overall	0.862	0.931	0.961	1.904	1.958

Table 4 .

 4 [START_REF] Torczon | On the convergence of pattern search algorithms[END_REF] Success rate (SR) and speedup performance (SPU) comparison of KA-CMA-ES using Fixed Individual-based Control (FIC) and Mixed Individual-based Control (MIC).

	Function D λ	FIC-Mean SR SPU SR SPU SR SPU SR SPU SR SPU SR FIC-SD FIC-SLB FIC-POI FIC-EI MIC-Mean SPU	MIC-SLB SR SPU
		6	1	1.94	1	1.82	1	1.84	1	1.87	1	1.90	1	2.53	1	2.53
	f1	8 10 10	1 1	1.82 1.62	1 1	1.63 1.25	1 1	1.85 1.64	1 1	1.81 1.62	1 1	1.87 1.64	1 1	2.25 1.77	1 1	1.98 1.56
		20 12	1	1.34	1	1.00	1	1.32	1	1.31	1	1.32	1	1.38	1	1.27
		6	1	1.91	1	1.85	1	1.90	1	1.92	1	1.90	1	2.61	1	2.63
	f2	8 10 10	1 1	1.95 2.04	1 1	1.69 1.86	1 1	1.74 1.77	1 1	1.98 1.89	1 1	1.82 1.91	1 1	2.38 2.43	1 1	2.46 2.34
		20 12	1	1.91	1	1.81	1	1.95	1	1.97	1	2.04	1	2.30	1	2.04
		6	1	1.89	1	1.83	1	1.91	1	1.85	1	1.90	1	2.49	1	2.49
	f3	8 10 10	1 1	1.83 1.55	1 1	1.62 1.19	1 1	1.85 1.62	1 1	1.80 1.61	1 1	1.80 1.62	1 1	2.29 1.57	1 1	1.94 1.54
		20 12	1	1.14	1	1.04	1	1.22	1	1.24	1	1.25	1	1.04	1	1.19
		6	1	1.88	1	1.84	1	1.85	1	1.85	1	1.94	1	2.53	1	2.56
	f4	8 10 10	1 1	1.85 1.51	1 1	1.65 1.28	1 1	1.87 1.58	1 1	1.88 1.61	1 1	1.85 1.58	1 1	2.29 1.52	1 1	1.98 1.44
		20 12	1	1.19	1	1.08	1	1.27	1	1.30	1	1.28	1	1.06	1	1.19
		6	1	1.76	1	1.76	1	1.79	1	1.82	1	1.80	1	2.43	1	2.41
	f5	8 10 10	1 1	1.81 1.88 0.92 1.55 0.96 1.80 0.96 1.75 0.96 1.82 0.84 1.88 0.84 1.79 1 1.71 1 1.84 1 1.87 1 1.77 1 2.30 1 2.27
		20 12 0.72 1.79 0.8 1.67 0.6 1.51 0.76 1.86 0.84 1.99 0.6	1.73	0.6 1.72
		6	1	1.91	1	1.90	1	1.95	1	1.84	1	1.95 0.96 2.41	1	2.29
	f6	8 0.96 1.24	1	1.30	1	1.67	1	1.77	1	1.68	1	1.24	1	1.77
		10 10 0.16 0.14 0.92 0.68 0.92 0.96 0.8 0.82	1	1.13	0	0.0	0.96 1.22
		6	1	1.84	1	1.98	1	1.91	1	1.82	1	1.94	1	2.50	1	2.37
	f7	8 0.96 1.92 0.96 1.76 0.92 1.84 0.96 1.91 0.92 1.80 0.92 2.17 0.92 1.97
		10 10 0.84 1.51 0.92 1.67 0.88 1.78 0.76 1.63 0.96 1.95 0.84 1.36 0.96 1.92
		6	1	1.87	1	1.88	1	1.92 0.92 1.79	1	1.92 0.88 2.29 0.96 2.33
	f8	8 0.96 1.99 0.96 1.71 0.96 1.97	1	2.09	1	2.06 0.92 2.44	1	2.42
		10 10 0.84 1.43 0.92 1.15 0.84 1.45 0.76 1.31	1	1.69 0.56 1.05 0.96 1.71
		12	1	1.94	1	1.85	1	1.97	1	1.86	1	1.90	1	2.23	1	2.08
	f9	16	1	2.00	1	1.83	1	2.02 0.96 1.89	1	1.99	1	2.54 0.96 2.23
		10 20 0.92 1.72 0.88 1.27 0.88 1.63 0.92 1.76 0.88 1.65 0.88 1.92 0.92 1.84
		12 0.96 1.72	1	1.82	1	1.89	1	1.92	1	1.88 0.64 1.38	1	2.10
	f10	16 0.96 1.90	1	1.70	1	2.03	1	1.95	1	1.97 0.84 2.14 0.96 2.33
		10 20 0.64 1.28 0.92 1.27 0.68 1.36 0.88 1.74 0.92 1.78 0.2	0.46 0.52 1.17
		12 0.96 1.87	1	1.91 0.96 1.90	1	1.98	1	2.03 0.96 2.21	1	2.27
	f11	16 0.84 1.70 0.76 1.42	1	2.06 0.88 1.76 0.8 1.63 0.84 2.21	0.8 2.02
		10 20 0.56 1.35 0.8 1.44 0.84 2.07 0.56 1.37 0.96 2.37 0.52 1.47 0.64 1.72
	f12	50 0.8 1.74 0.96 1.91 0.92 2.00 0.84 1.76 0.92 1.97 0.84 2.01 0.84 2.04 140 0.4 1.05 0.56 1.02 0.4 0.94 0.72 1.68 0.72 1.65 0.24 0.82 0.4 1.17

Table 4 .

 4 [START_REF] Audet | Mesh Adaptive Direct Search Algorithms for Constrained Optimization[END_REF] Average success rate and speedup performance of KA-CMA-ES using Fixed and Mixed Individual-based Control (FIC and MIC)

	Category	Average Success Rate (SR) FIC-Mean FIC-SD FIC-SLB FIC-POI FIC-EI MIC-Mean MIC-SLB
	D=2	0.977	0.997	0.990	0.980	0.993	0.940	0.983
	D=5	0.923	0.937	0.940	0.960	0.953	0.897	0.920
	D=10	0.815	0.935	0.909	0.876	0.971	0.713	0.891
	D=20	0.944	0.960	0.920	0.952	0.968	0.920	0.920
	Unimodal	0.948	0.982	0.972	0.971	0.988	0.929	0.972
	Multimodal	0.846	0.911	0.891	0.889	0.943	0.737	0.854
	Overall	0.912	0.957	0.944	0.942	0.972	0.862	0.931
	Category	Average Speedup Performance (SPU) FIC-Mean FIC-SD FIC-SLB FIC-POI FIC-EI MIC-Mean MIC-SLB
	D=2	1.857	1.862	1.902	1.858	1.920	2.302	2.341
	D=5	1.757	1.587	1.806	1.866	1.824	2.090	2.046
	D=10	1.456	1.327	1.606	1.554	1.738	1.403	1.661
	D=20	1.473	1.319	1.455	1.536	1.576	1.503	1.483
	Unimodal	1.661	1.554	1.701	1.716	1.747	1.941	1.957
	Multimodal	1.683	1.584	1.800	1.775	1.892	1.798	1.960
	Overall	1.669	1.564	1.736	1.737	1.798	1.891	1.958

Table 4 .

 4 [START_REF] Jones | Lipschitzian optimization without the Lipschitz constant[END_REF] Results of KA-CMA-ES using Approximate Ranking Procedure (ARP)

	Function	D	λ	SR	ARP-Mean SP	SPU	SR	ARP-EI SP	SPU	CMA-ES SR SP
		2	6	1	82	4.04	1	86	3.85	1	331
	f1	5 10	8 10	1 1	351 1053	2.53 1.69	1 1	352 1014	2.52 1.76	1 1	888
		20	12	1	2453	1.34	1	2444	1.35	1
		2	6	1	84	4.61	1	87	4.45	1	387
	f2	5 10	8 10	1 1	345 1067	3.27 2.76	1 1	365 996	3.09 2.95	1 1
		20	12	1	2530	2.95	1	2563	2.91	1
		2	6	1	84	4.24	1	90	3.96	1	356
	f3	5 10	8 10	1 1	393 1272	2.49 1.63	1 1	387 1233	2.53 1.68	1 1	980
		20	12	1	3875	1.12	1	3540	1.22	1
		2	6	1	100	4.17	1	104	4.01	1	417
	f4	5 10	8 10	1 1	490 1832	2.50 1.56	1 1	484 1667	2.53 1.71	1 1
		20	12	1	6660	1.12	1	5815	1.28	1
		2	6	1	76	3.64	1	74	3.74	1	277
	f5	5 10	8 10	1 0.96	424 1452	2.71 2.45	1 0.96	397 1471	2.89 2.42	1 1
		20	12	0.44	9849	1.42	0.64	6437	2.18	0.76	14014
		2	6	1	258	3.02	1	193	4.03	1	778
	f6	5	8	0.8	1462	1.54	1	1225	1.84	1
		10	10	0.12	41047	0.12	0.88	6148	0.83	1
		2	6	1	172	3.85	1	157	4.22	1	663
	f7	5	8	0.84	1163	2.12	0.92	968	2.55	0.96
		10	10	0.84	5658	1.46	0.88	4249	1.94	0.8
		2	6	0.84	225	3.19	0.96	191	3.75	1	717
	f8	5	8	0.8	906	2.25	1	726	2.81	0.92
		10	10	0.88	2400	1.56	0.96	2164	1.73	0.96
		2	12	1	150	3.24	1	149	3.26	1	486
	f9	5	16	0.88	558	2.67	0.96	501	2.97	0.96
		10	20	0.92	1717	1.78	1	1585	1.92	1
		2	12	0.4	980	1.08	1	368	2.86	1
	f10	5	16	0.64	1813	1.73	0.96	1146	2.74	0.96
		10	20	0.56	6724	1.03	0.96	3903	1.77	0.92
		2	12	0.88	199	3.00	0.96	182	3.28	1	597
	f11	5	16	0.72	752	2.26	0.96	569	2.99	0.96
		10	20	0.8	2264	1.97	0.8	2260	1.97	0.76
	f12	2 5	50 140	0.88 0.72	1042 6363	2.28 1.57	0.96 0.8	916 5884	2.59 1.70	0.68 0.8	10020

Table 4 .

 4 [START_REF] Huyer | Global Optimization by Multilevel Coordinate Search[END_REF] Average success rate and speedup of KA-CMA-ES using ARP. For generation-based control, the Fixed Generation-based Control (FGC) where one generation is controlled in each two generations and the proposed Adaptive Generationbased Control (AGC) are investigated. The results of generation-based control are listed in

	Category	Average Success Rate (SR) ARP-Mean ARP-EI CMA-ES ARP-Mean Average Speedup (SPU) ARP-EI
	D=2	0.917	0.990	0.973	3.362	3.667
	D=5	0.867	0.967	0.963	2.304	2.597
	D=10	0.825	0.949	0.949	1.637	1.881
	D=20	0.888	0.928	0.952	1.590	1.788
	Unimodal	0.923	0.972	0.982	2.475	2.633
	Multimodal	0.780	0.949	0.923	2.114	2.597
	Overall	0.873	0.964	0.961	2.349	2.620
	4.5.5 Experiments on KA-CMA-ES using Generation-based Control

Table 4 .

 4 17. The average success rate and speedup performance are presented in Table 4.18 and Figure 4.10. In these tables, FGC stands for fixed generation-based control, AGCw Q

	indicates the adaptive generation-based control based on model quality w Q ,	selection AGC-Q
	signifies adaptive generation-based control based on model quality measure selection Q	, and
	rank AGC- represents AGC based on model quality measure rank . 
	In Table 4.17, we mark the SPU that do not exceed one ( SPU 1  ) in red bold format,
	which means the corresponding algorithm does not improve the success performance on the
	corresponding test problem. In other words, we can say that the investigated algorithm
	performs worsen than standard CMA-ES on the test problem if the SPU of the algorithm on
	a test problem is not larger than 1. For FGC, the number of problems on which SPU is not
	larger than one is 6. This number for AGC-w Q is 5, for	selection AGC-Q	is 7, and for	AGC-	rank
	is 5. The smaller this number is, the corresponding algorithm is more stable and reliable.
	From this perspective, AGC-w Q and	rank AGC- are slightly more stable and reliable than
	FGC and	selection AGC-Q	.

Table 4 .

 4 [START_REF] Jin | Fitness Approximation in Evolutionary Computation: A Survey[END_REF] Results of KA-CMA-ES using FGC and AGC

	Function D	λ	SR	FGC SP	AGC-w Q SPU SR SP SPU	selection AGC-Q SR SP SPU SR	rank SP SPU SR AGC- CMA-ES SP
		2	6	1	179 1.85 1 148	2.24	1	151 2.19	1	156 2.12 1	331
	f1	5 10	8 10	1 1	523 1.70 1 499 1308 1.36 1 1264	1.78 1.41	1 1	473 1.88 1289 1.38	1 1	586 1.52 1 1589 1.12 1	888
		20	12	1	3127 1.05 1 3205	1.03	1	3202 1.03	1	3376 0.98 1
		2	6	1	213 1.82 1 164	2.36	1	164 2.36	1	170 2.28 1	387
	f2	5 10	8 10	1 1	654 1.72 1 477 1675 1.76 1 1181	2.36 2.49	1 1	460 2.45 1168 2.52	1 1	523 2.16 1 1523 1.93 1
		20	12	1	4584 1.63 1 4068	1.84	1	4018 1.86	1	4820 1.55 1
		2	6	1	195 1.83 1 158	2.25	1	153 2.33	1	172 2.07 1	356
	f3	5 10	8 10	1 1	588 1.67 1 549 1635 1.27 1 1630	1.79 1.27	1 1	533 1.84 1616 1.28	1 1	650 1.51 1 1925 1.08 1	980
		20	12	1	4958 0.87 1 4914	0.88	1	4838 0.89	1	4523 0.96 1
		2	6	1	232 1.80 1 179	2.33	1	176 2.37	1	195 2.14 1	417
	f4	5 10	8 10	1 1	724 1.69 1 672 2295 1.25 1 2361	1.82 1.21	1 1	658 1.86 2384 1.20	1 1	793 1.54 1 2636 1.08 1
		20	12	1	8347 0.89 1 8530	0.87	1	8030 0.93	1	7536 0.99 1
		2	6	1	155 1.79 1 134	2.07	1	128 2.16	1	145 1.91 1	277
	f5	5 10	8 10 0.8 2588 1.37 0.92 1844 1 667 1.72 0.96 584	1.97 1.93	1 0.88 1978 1.80 0.92 2530 1.41 1 561 2.05 1 659 1.74 1
		20	12 0.64 10132 1.38 0.72 7362	1.90	0.64 8377 1.67 0.8 9835 1.42 0.76 14014
		2	6	1	444 1.75 1 349	2.23	1	353 2.20	1	432 1.80 1	778
	f6	5	8	1	2098 1.07 0.84 2342	0.96	0.84 2182 1.03	1	1802 1.25 1
		10	10 0.04 1.19E+05 0.04 0		0.0	0	 0.0 1 4881 1.04 1
		2	6	1	375 1.77 1 296	2.24	1	296 2.24	1	309 2.15 1	663
	f7	5	8 0.88 1594 1.55 0.84 1642	1.50	0.8 1692 1.46 0.88 1762 1.40 0.96
		10	10 0.92 6283 1.31 0.88 6980	1.18	0.96 6248 1.32 0.92 6686 1.23 0.8
		2	6	1	401 1.79 0.96 329	2.18	0.96 337 2.13	1	377 1.90 1	717
	f8	5	8	1	1091 1.87 0.96 1045	1.95	0.68 1451 1.41 0.96 1288 1.58 0.92
		10	10 0.64 4173 0.90 0.8 3422	1.10	0.64 4314 0.87	1	3181 1.18 0.96
		2	12	1	256 1.90 1 214	2.27	1	216 2.25	1	253 1.92 1	486
	f9	5	16 0.88 896 1.66 1 643	2.32	0.96 685 2.17 0.96 1058 1.41 0.96
		10	20 0.8 2241 1.36 0.8 2094	1.46	0.8 2115 1.44 0.96 2880 1.06 1
		2	12 0.92 665 1.58 0.84 631	1.67	0.96 555 1.90 0.96 908 1.16 1
	f10	5	16 0.8 2207 1.42 0.48 3054	1.03	0.4 3841 0.82	1	2714 1.16 0.96
		10	20 0.24 16028 0.43 0		0.0	0.2 18420 0.37 0.88 6927 1.00 0.92
		2	12 0.96 322 1.85 0.88 275	2.17	0.96 256 2.33	1	299 2.00 1	597
	f11	5	16 0.72 1208 1.41 0.76 926	1.84	0.72 997 1.71 0.88 1271 1.34 0.96
		10	20 0.64 3146 1.42 0.52 3530	1.26	0.4 4690 0.95 0.64 4814 0.93 0.76
	f12	2 5	50 0.92 1186 2.00 0.8 1184 140 0.28 14143 0.71 0.12 27611 0.36 2.00	0.92 1033 2.30 0.84 1219 1.95 0.68 0.44 7839 1.28 0.8 7858 1.28 0.8 10020

Table 4 .

 4 [START_REF] Jin | A comprehensive survey of fitness approximation in evolutionary computation[END_REF] Average success rate and speedup of KA-CMA-ES using FGC and AGC. Figure 4.10 Average success rate and speedup of KA-CMA-ES using FGC and AGC.) has the highest average success rate on all the problems. On multimodal problems, the average success rate of FGC, AGC-

	Category	Average Success Rate (SR) FGC AGC-w Q selection AGC-Q rank AGC-	Average Speedup (SPU) CMA-ES FGC AGC-w Q selection AGC-Q AGC-	rank
	D=2	0.983 0.957	0.983	0.983	0.973 1.810 2.168	2.230	1.949
	D=5	0.880 0.830	0.820	0.957	0.963 1.516 1.640	1.662	1.490
	D=10 0.735 0.720	0.716	0.938	0.949 1.133 1.210	1.194	1.187
	D=20 0.928 0.944	0.928	0.960	0.952 1.166 1.304	1.276	1.179
	Unimodal 0.934 0.929	0.928	0.982	0.982 1.458 1.689	1.704	1.553
	Multimodal 0.771 0.709	0.717	0.920	0.923 1.450 1.543	1.566	1.418
	Overall 0.877 0.852	0.854	0.960	0.961 1.455 1.638	1.656	1.505

From the radar charts in

Figure 4.10, rank AGC- (adaptive generation control using rank 

Table 4 .

 4 [START_REF] Jin | Surrogate-assisted evolutionary computation: Recent advances and future challenges[END_REF] Performance comparison of KA-CMA-ES algorithms, including PS-Interval, CPSw Q ,

			rank CPS- , FIC-EI, MIC-SLB, ARP-EI, AGC-Q and w	rank AGC- .
	Function D λ	PS-Interval CPS-Q w SR SPU SR SPU SR SPU SR SPU SR SPU SR SPU SR SPU SR rank CPS- FIC-EI MIC-SLB ARP-EI AGC-w Q AGC- SPU rank
		2 6	1	1.29 1 1.61 1	1.64 1 1.90 1 2.53 1 3.85 1 2.19	1	2.12
	f1	5 8 10 10	1 1	1.26 1 1.64 1 1.27 1 1.59 1	1.63 1 1.87 1 1.98 1 2.52 1 1.88 1.46 1 1.64 1 1.56 1 1.76 1 1.38	1 1	1.52 1.12
		20 12	1	1.13 1 1.19 1	1.15 1 1.32 1 1.27 1 1.35 1 1.03	1	0.98
		2 6	1	1.21 1 1.51 1	1.50 1 1.90 1 2.63 1 4.45 1 2.36	1	2.28
	f2	5 8 10 10	1 1	1.12 1 1.42 1 1.17 1 1.58 1	1.39 1 1.82 1 2.46 1 3.09 1 2.45 1.65 1 1.91 1 2.34 1 2.95 1 2.52	1 1	2.16 1.93
		20 12 0.96 1.13 1 1.98 1	1.73 1 2.04 1 2.04 1 2.91 1 1.86	1	1.55
		2 6	1	1.30 1 1.62 1	1.60 1 1.90 1 2.49 1 3.96 1 2.33	1	2.07
	f3	5 8 10 10	1 1	1.29 1 1.66 1 1.25 1 1.54 1	1.65 1 1.80 1 1.94 1 2.53 1 1.84 1.41 1 1.62 1 1.54 1 1.68 1 1.28	1 1	1.51 1.08
		20 12	1	1.03 1 0.82 1	1.05 1 1.25 1 1.19 1 1.22 1 0.89	1	0.96
		2 6	1	1.28 1 1.64 1	1.65 1 1.94 1 2.56 1 4.01 1 2.37	1	2.14
	f4	5 8 10 10	1 1	1.30 1 1.72 1 1.24 1 1.36 1	1.69 1 1.85 1 1.98 1 2.53 1 1.86 1.29 1 1.58 1 1.44 1 1.71 1 1.20	1 1	1.54 1.08
		20 12	1	1.09 0.6 0.49 1	1.10 1 1.28 1 1.19 1 1.28 1 0.93	1	0.99
		2 6	1	1.28 1 1.70 1	1.70 1 1.80 1 2.41 1 3.74 1 2.16	1	1.91
	f5	5 8 10 10	1 1	1.26 1 2.01 1 1.27 1 2.24 1	1.93 1 1.77 1 2.27 1 2.89 1 2.05 2.08 0.96 1.82 0.84 1.79 0.96 2.42 0.88 1.80 0.92 1.41 1 1.74
		20 12 0.84 1.34 1 3.28 1	2.89 0.84 1.99 0.6 1.72 0.64 2.18 0.64 1.67 0.8	1.42
		2 6	1	1.35 1 1.74 1	1.67 1 1.95 1 2.29 1 4.03 1 2.20	1	1.80
	f6	5 8	1	1.41 1 1.82 1	1.66 1 1.68 1 1.77 1 1.84 0.84 1.03	1	1.25
		10 10	1	1.37 1 1.76 1	1.54 1 1.13 0.96 1.22 0.88 0.83 0	0.0	1	1.04
		2 6	1	1.48 1 2.15 1	2.07 1 1.94 1 2.37 1 4.22 1 2.24	1	2.15
	f7	5 8 0.96 1.47 0.96 2.24 0.96 2.02 0.92 1.80 0.92 1.97 0.92 2.55 0.8 1.46 0.88 1.40
		10 10 0.96 1.63 0.96 2.06 0.96 1.81 0.96 1.95 0.96 1.92 0.88 1.94 0.96 1.32 0.92 1.23
		2 6	1	1.29 1 1.66 0.96 1.55 1 1.92 0.96 2.33 0.96 3.75 0.96 2.13	1	1.90
	f8	5 8	1	1.44 1 1.96 1	1.91 1 2.06 1 2.42 1 2.81 0.68 1.41 0.96 1.58
		10 10	1	1.34 1 1.79 1	1.59 1 1.69 0.96 1.71 0.96 1.73 0.64 0.87	1	1.18
		2 12	1	1.29 1 1.42 1	1.41 1 1.90 1 2.08 1 3.26 1 2.25	1	1.92
	f9	5 16	1	1.31 1 1.52 1	1.51 1 1.99 0.96 2.23 0.96 2.97 0.96 2.17 0.96 1.41
		10 20	1	1.25 1 1.49 1	1.39 0.88 1.65 0.92 1.84 1 1.92 0.8 1.44 0.96 1.06
		2 12 0.96 1.17 0.52 0.71 0.92 1.22 1 1.88 1 2.10 1 2.86 0.96 1.90 0.96 1.16
	f10	5 16	1	1.33 1 1.65 1	1.47 1 1.97 0.96 2.33 0.96 2.74 0.4 0.82	1	1.16
		10 20	1	1.40 1 1.75 1	1.46 0.92 1.78 0.52 1.17 0.96 1.77 0.2 0.37 0.88 1.00
		2 12	1	1.30 1 1.50 1	1.44 1 2.03 1 2.27 0.96 3.28 0.96 2.33	1	2.00
	f11	5 16	1	1.31 0.96 1.50 0.92 1.41 0.8 1.63 0.8 2.02 0.96 2.99 0.72 1.71 0.88 1.34
		10 20 0.92 1.51 0.92 1.81 0.8 1.50 0.96 2.37 0.64 1.72 0.8 1.97 0.4 0.95 0.64 0.93
	f12	2 50 0.92 1.59 0.8 1.40 0.96 1.67 0.92 1.97 0.84 2.04 0.96 2.59 0.92 2.30 0.84 1.95 5 140 0.84 1.34 0.64 1.10 0.84 1.39 0.72 1.65 0.4 1.17 0.8 1.70 0.44 1.28 0.8 1.28

Table 4 .

 4 [START_REF] Shi | A Survey of Fitness Approximation Methods Applied in Evolutionary Algorithms[END_REF] Comparison of Average success rate (SR) and speedup performance (SPU) of KA-CMA-ES algorithms, including PS-Interval, CPS-w Q ,

	rank CPS- , FIC-EI, MIC-SLB, ARP-EI, AGC-w Q and
	rank AGC- .

Figure 4.11 Comparison of average success rate (SR) of KA-CMA-ES algorithms.

  Generally, optimization algorithms can be divided into two main categories: nongradient methods (or called direct search methods), which do not employ gradients and involve only evaluations of the objective function, and gradient-based methods, which require computations of the derivatives of the objective function[99]. Gradient-based methods, generally, require a much smaller number of design cycles to converge to an optimum compared to non-gradient methods. However, only convergence to a local minimum is guaranteed for gradient-based methods, while non-gradient methods are able to find global minimum. In addition, gradient-based methods are limited in cases that the objective function is always differentiable and has continuous derivatives over the design domain due to the requirements of computations of derivatives of the objective function.

Table 5 .

 5

	1 Equations of elastic-plastic damage model
	Strain tensor additive split

  By introducing Equation (5.23) and (5.20) 1 into the definition of the damage energy release rate of Equation (5.11), the Equation (5.18) 3 can be expressed as:

	   	,	 trial 11 n Dq n  		11 11 3 p y n nn G DD           		0.	(5.24)
	1 nn DD   	1	 11 , nn YD Dr 1 n           	s		0,	(5.25)
	where
	 11 ,.      2 2 trial 11 1 2 13 nn n 6 (1 ) nn 2 1 n

  and (5.25), which is written as

	trial 11 n	3 11 p y n 11 nn		0,
	21 nn	1	n	1	 11 , nn		s	0.

Table 5 .

 5 2 Stress update algorithm for elastic-plastic damage model

	(i) Elastic predictor
		e n	trial 11 trial ; e p n n        	p n
	s	trial 1 n 	trial d 1 2; ee trial 1 v 1 trial n n n G p K     
	trial 1 n q 		3 2	trial 1 n  ss :	trial 1 n 
	(ii) Check plastic admissibility
	If	trial  	trial 1 n q 		  p n  y	0  Then
			Set   	  11 trial nn   	and Exit
	Else go to (iii)
	(iii) Return-mapping (solve the system for the unknowns   and	1 D  ) n
			trial 11 n		3 11 p y n 11 nn	0
			21 nn	1	n	1	 11 , nn		s	0

Table 5 .

 5 [START_REF] Bä Ck | Contemporary Evolution Strategies[END_REF] Results and computational costs of strain hardening parameter identification

	Parameter	A		0	n	FES
	Lower Bound	600	0	0.1	-
	Upper Bound	900	0.03	0.4	-
	The Standard CMA-ES	754.1402	0.0087	0.2264	533
	KA-CMA-ES using ARP-EI	754.2124	0.0087	0.2264	211

Table 5 .

 5 [START_REF] Koziel | Computational Optimization[END_REF] Results and computational costs of damage parameter identification

	5. Applications in Material Parameter Identification

Résumé de la Thèse en Français 6.1 Introduction Générale

  This chapter presents a ductile damage model and applies a previously proposed KA-CMA-ES algorithm (KA-CMA-ES using ARP-EI) in inverse method of parameter identification. The ductile damage model, which combines the von Mises yield criterion, Swift's hardening law and Lemaitre's damage model, is implemented by the fully implicit elastic predictor/return-mapping scheme in ABAQUS through the subroutine VUMAT. The numerical implementation algorithm is comprehensively detailed in this chapter. In order to improve the efficiency of the inverse method for parameter identification, previously proposed KA-CMA-ES (ARP-EI) algorithm is used. The proposed ductile damage model is employed for A 2017-T4 and the inverse method using ARP-EI is applied to identify the material parameters based on standard tensile test data.En raison de la forme implicite de la fonction objectif, les algorithmes é volutifs (EAs), qui sont une classe d'optimiseurs globaux sans dé rivé es et puissants, sont approprié s pour ré soudre des problè mes d'optimisation coûteux. Cependant, la principale difficulté à utiliser les EAs est que ces derniers ont habituellement besoin d'un grand nombre d'é valuation de la fonction objectif avant d'obtenir un ré sultat satisfaisant. Par consé quent, les EAs assisté es par modè le de substitution ont é té motivé s par la ré duction des coûts du temps de calcul.

			Min	f	  x
			s.c.	 l x u
	où   : d f x		est la fonction objectif qui est é valué e en exé cutant la simulation
	L'optimisation est largement demandé e et appliquée en science de l'ingénieur. Motivé es par les demandes industrielles et de la recherche scientifique, beaucoup de techniques d'optimisation ont é té dé veloppé es. L'algorithme choisi pour un problè me d'optimisation dé pend en grande partie du type de problè me, de la qualité de la solution souhaité e, des ressources informatiques disponibles, du temps CPU et de l'expertise des dé cideurs.  Pour de nombreux problè mes d'optimisation de CAO, la nature de la fonction objectif peut ê tre non lisse, multimodale, discontinue et mal conditionné e. Ces difficulté s entraî nent un coût de calcul é levé pour trouver l'optimum. Nous nous concentrons sur les problè mes d'optimisation coûteux sous contraintes dé crits comme suit : numé rique, le vecteur   T 1 ,, d  x repré sente les d variables de conception, les vecteurs l et u sont les limites infé rieure et supé rieure, respectivement. Les contraintes limitent xx l'espace de recherche à       1 1 2 2 , , , dd l u l u l u    .

The results show that Swift's law is adequate to describe the strain hardening behavior of A 2017-T4. With the incorporation of Lemaitre's ductile damage into the elastoplastic model, ductile damage effect can be modeled. The mechanical behavior of A 2017-T4 under tension is accurately model by the presented elastic-plastic damage model.

Applications of the KA-CMA-ES algorithm prove that this algorithm is useful and promising for parameter identification. It enhances the efficiency of the parameter identification process and also gives reliable results.

6.

6.1.1 Motivation

Une problè me d'optimisation est dit difficile s'il est coûteux en temps CPU. Ce coût est gé né ralement dû:  La fonction objectif est é valué e sur la base d'une simulation numé rique, qui peut ê tre coûteuse (requiert de quelques minutes à des heures voire des jours de temps de calcul pour chaque cycle de simulation).

 Il n'existe pas d'expression analytique explicite pour la fonction objectif ou ses dé rivé es. Ainsi, les algorithmes sans dé rivation sont né cessaires pour ré soudre ce type de problè me d'optimisation. Ces derniers né cessitent plusieurs é valuations de la fonction objectif que celles utilisant les dé rivé s.

6.1.2 Objectif

L'objectif de cette thè se est de dé velopper des algorithmes d'optimisation puissants qui peuvent traiter plus efficacement les problè mes d'optimisation coûteux. Ce travail se concentre sur la straté gie d'é volution (ES) assisté e par modè le de substitution pour des problè mes d'optimisation coûteux. La straté gie d'évolution avec adaptation de matrice de covariance (CMA-ES) et le modè le de krigeage sont choisis comme les deux composantes de la straté gie d'é volution assisté e par modè le de substitution. Cette CMA-ES assisté e par le modè le de krigeage est abré gé e par KA-CMA-ES dans ce manuscrit. Notre objectif est d'é tudier la CMA-ES existante et de dé velopper de nouveaux algorithmes efficaces de KA-CMA-ES pour des problè mes coûteux.

6.2.2 Deuxième Chapitre: Les Stratégies d'Evolution

	Ce chapitre donne une description complè te des straté gies d'é volution (ES). Une brè ve
	introduction des ES est fournie en premier lieu. Ensuite, les principes et les opé rateurs
	é volutifs utilisé s dans les ES, c'est-à -dire la sé lection, la recombinaison et la mutation, sont
	dé crits dans ce chapitre. Par la suite, nous pré sentons trois algorithmes typiques de ES,
	(1+1)-ES, 	 , -ES    I	avec adaptation cumulative de taille d'é chelon et straté gie
	d'é volution avec l'adaptation de matrice de covariance (CMA-ES).
	Algorithme 6.1 Le  1: Donné es: dd , -CMA-ES     w   , 4 3ln	,	2

would express my sincere acknowledgment to China Scholarship Council (CSC) for its financial support of my study in France and to my supervisor Professor Abdelkhalak EL HAMI for his guidance, patience, encouragement and providing me grateful for the period of time when we worked together and discussed about our research papers. I would also like to express my sincere gratitude to all the members of my doctoral thesis committee for their insightful comments on my thesis and encouragement.

To all my colleagues in LMN and my friends who study in INSA de Rouen and CORIA,

Kriging-Assisted CMA Evolution Strategy

This chapter focuses on Kriging-assisted evolution strategies, under the framework of surrogate-assisted evolutionary optimization. A brief introduction of surrogate-assisted evolution strategy is presented at Section 4.

Introduction

For continuous optimization problems, evolution strategy (ES) generally works better than other evolutionary algorithms, such as genetic algorithm (GA) [67]. In recent black-box optimization benchmarking (2009 and 2010 GECCO), CMA-ES shows its outstanding performance and has proven to be one of the best-performing search strategies for realvalued black-box optimization [START_REF] Du | Evolutionary Strategies[END_REF]. Therefore, in this thesis, we adopt ES to solve expensive optimization problems. However, ES, like other population based EAs, require a large number of fitness function evaluations before obtaining a satisfying solution. Additionally, in expensive optimization problems of real world applications, such as engineering design optimization, no analytical fitness function exists and the evaluations of fitness function are by means of expensive numerical simulations or experiments. Thus, large number of fitness evaluations is not practical or prohibitive. Computational cost, consequently, has been a crucial challenge in application of ES (and other EAs) to expensive optimization problems. performance. Particularly, MIC-SLB has the highest average speedup performance and acceptable success rate. Furthermore, in KA-CMA-ES using the modified approximate ranking procedure (ARP), the Mean and EI metrics have been studied. The results show that ARP-EI clearly outperforms ARP-Mean. This also proves that EI is the most promising metric in evolution control. Therefore, it can be stated that the ARP-EI is outstanding among all the investigated KA-CMA-ES algorithms, according to the speedup performance or success performance.

The results of experiments of

Applications in Material Parameter

Identification

This chapter presents applications of the proposed KA-CMA-ES algorithm in material parameter identification, which involves expensive optimization problems. An elasticplastic damage model, in which material's strain hardening behavior is described by Swift law and ductile damage is modeled by Lemaitre's damage model, is presented and implemented in numerical simulation through VUMAT subroutine in ABAQUS. Then, the KA-CMA-ES using Approximate Ranking Procedure with metric EI (ARP-EI) is applied in inverse method of material parameter identification.

Introduction

Constitutive model and parameter identification can be seen as the two main aspects of material modeling. Specifically, constitutive model describes material's behavior by the use of mathematical equations or formulations under the framework of constitutive theories; while parameter identification determines the unknown parameters in the constitutive model on the basis of experimental data. Both the constitutive model and parameter identification play important roles in material modeling, and neither of the two aspects could be neglected.

Constitutive Model

The material's constitutive model has experienced huge development during last decades with the improvements both in constitutive theories and their applications. The descriptions of elasticity and plasticity have been comprehensively investigated and theorized in elastic-plastic theories (or elastoplastic theories). Even the more complicated material's behavior, such as viscoplasticity, fatigue, damage and anisotropy, can be modeled at present. With regard to damage, which is a critical and popular issue in engineering of the void development are derived by means of micromechanics analysis [94]. The Gurson-Tvergaard-Needleman (GTN) model is the representative model of this approach. The second approach to ductile damage analysis, called phenomenological approach, is based on the CDM theory proposed by Lemaitre [START_REF] Lemaitre | A three-dimensional ductile damage model applied to deep-drawing forming limits[END_REF], who has described the ductile damage phenomenon by a phenomenological model with the combination of continuum mechanics and irreversible thermodynamics. Using the method of local state and internal variables, Lemaitre's model strongly coupled damage and elastoplasticity at the constitutive level. In this paper, Lemaitre's model is used to study the ductile damage of the material.

Material Parameter Identification

In the aspect of parameter identification, inverse techniques have become popular in nowadays. In the past, the hand fitting method and the trial and error method are commonly used in material parameter identification [START_REF] Stein | Encyclopedia of computational mechanics[END_REF]. However, the inverse method has become a promising method since optimization techniques and simulation are widely applied. The inverse method considers parameter identification as an optimization problem and generally gives satisfy results [97]. This method tries to find a set of parameters which yields the simulated responses as closely as possible to the experimental response, in other words, inverse method aims at obtaining the set of material parameters by solving the optimization problem which minimize the difference between simulation results and experiment data.

Commonly, in the process of parameter identification by inverse method, the parameters are driven by optimization algorithm, and the objective function, which measures the difference between the simulated and experimental responses, is repeatedly evaluated by numerical simulation until the computational scheme converges. Since optimization algorithm (or optimization technique), which is used to find the solution to the optimization problem, has influences on the convergence speed and the solution of the optimization problem [START_REF] Hami | Uncertainty and Optimization in Structural Mechanics[END_REF], appropriate algorithm should be chosen in the inverse problem of parameter identification.

Conclusions and Perspectives

In this thesis, we have comprehensively explored the Kriging-Assisted Covariance Matrix Adaptation Evolution Strategy (KA-CMA-ES) for expensive optimization problems.

Furthermore, engineering application of the developed new algorithm is performed in material parameter identification. This final chapter gives a summary of what we have learned from our works and points out the most promising direction for future research.

Conclusions

We Additionally, the performance of KA-CMA-ES using PS and CPS are compared. It is apparent that pre-selection with model impact control outperforms pre-seleciton without model impact control for most of the test problems according to the speedup performance. Thus, model impact control is suggested in KA-CMA-ES using pre-selection.

For KA-CMA-ES using fixed individual-based control (FIC), five metrics (Mean, SD, SLB, POI, and EI) are studied. Among the five metric, these metric (SLB, POI and EI) which balances the exploitation of the surrogate and exploration of search space performs better than metrics that only consider exploitation or exploration (Mean, and SD). Take both SR and SP into consideration, FIC-EI has advantages over other metrics.

In our proposed mixed individual-based control (MIC), two possible metrics, Mean and SLB, are examined. The experimental results proved that MIC-SLB has higher success rate than that of MIC-Mean, especially on multimodal problems. Furthermore, generally, MIC-SLB has better average speedup performance than MIC-Mean. Taking both the success rate and speedup performance into consideration, it can be concluded that MIC-SLB outperforms MIC-Mean. From the comparison of KA-CMA-ES using FIC and MIC, the FIC-EI generally has the highest success rate. Other individual-based control algorithms have moderate 0. Conclusions and Perspectives 203 average success rate larger than 0.85. However, KA-CMA-ES using MIC significantly outperform FIC according to the speedup performance and success performance. Particularly, MIC-SLB has the highest average SPU and acceptable success rate. Therefore, MIC-SLB is preferable among algorithms of KA-CMA-ES using individual-based control.

We have modified the approximate ranking procedure and then embedded it into the CMA-ES, which is referred to as KA-CMA-ES using approximate ranking procedure (ARP).

Besides the Mean metric, the previously suggested metric EI is used in KA-CMA-ES using ARP. It has been concluded that ARP-EI outperforms ARP-Mean, considering both the success rate and success performance.

Moreover, the KA-CMA-ES using generation-based control has been investigated, where the fixed generation-based control (FGC) and the proposed adaptive generation-based control (AGC) are examined and compared. Generally, the KA-CMA-ES using AGC performs better than that using FGC. The adaptive generation control using rank The application of KA-CMA-ES using ARP-EI in inverse method of parameter identification for the elastic-plastic damage model also shows that the KA-CMA-ES improve the efficiency of CMA-ES significantly and thus greatly reduce the computational cost of material parameter identification.

Perspectives

The developed KA-CMA-ES algorithms have apparently improve the efficiency of the standard CMA-ES. However, there are still some open questions. Here we would like to mention several issues of future research, from our view of points.

Firstly, with respect to the surrogate models used in surrogate-assisted evolutionary computation, there is still no clear conclusion on which model is better than others. Thus, more surrogate models need to be investigated.

Regards to the KA-CMA-ES developed in this work, the experimental study and application on high dimensional problems (larger than 20 dimension) could be carried out in future works.

For high dimensional problems, the computational cost of surrogate model training would increase apparently. This may become a disadvantage of surrogate-assisted evolution algorithms. Therefore, surrogate-assisted evolutionary optimization for high dimensional problem is a topic worthy of studying.