Development of MARS - Multipoint Approximation method based on the Response Surface fitting

Motivation for the development of MARS: Design optimization and inverse problems with computationally expensive and noisy response functions.



Many real-life design optimization problems have the following common features:

Presence of numerical noise in the response functions means that function values, corresponding to two slightly different sets of design variables, may have a considerable random variation. Moreover, even function values, corresponding to the same set of design variables, may be different depending on the optimization search history and/or iteration history of an iterative response analysis. As an example, we mention FEM in conjunction with Adaptive Mesh Refinement (AMR), see Polynkine et al. (1996a), van Keulen et al. (1997) and van Keulen and Toropov (1997), for which the mesh density information can be influenced by the optimization history.

The level of noise can vary from negligible to significant (or even intolerable) depending on the formulation and the size of the optimization problem as well as the size and specific features of the response analysis model.

As an example of the noisy behavior of response functions, consider a simple optimization problem of a square plate with a circular central hole, Toropov et al. (1995). The radius of the hole is treated as a single design variable. The normalized stress constraint (Fig. 1) and normalized stability constraint (Fig. 2) are shown for a coarse and a fine mesh (Fig. 3 and 4). In each case the mesh density was kept constant. It can be seen that the the coarse mesh results possess a significantly higher level of noise. Another obvious observation is that the numerical solution is also characterized by an offset from the exact solution. The magnitude of the offset depends on the specified mesh density and tends to zero when the mesh density is increased.

Figure 1. Normalized stress constraint

Figure 2. Normalized stability constraint

Figure 3. Coarse mesh

Figure 4. Fine mesh

Another example (Polynkin et al. 1996b) relates to the optimization of a conical shell made of a continuous fibre reinforced thermoplastic material (Fig. 5). The response analysis in this case involves two modelling steps, namely simulation of the thermoforming process and the FE analysis of the obtained product. The problem is formulated as maximization of the strain energy subject to the constant volume constraint. The design variables are the height of the cone and the base radius. Fig. 6 shows the noisy behavior of the objective function.

Figure 5. CFRTP conical shell

Figure 6. Objective function for the CFRTP conical shell problem

Other design optimization problems exhibiting noisy behaviour of response functions have been considered by Free et al. (1987), Giunta et al. (1994), Narducci et al. (1995), Valorani and Dadone (1995), Toropov et al. (1996), Polynkine (1996). During the last decade, several approximation techniques were proposed which take account of the negative effect of noise.

Rasmussen (1990) investigated the influence of noise in the first order sensitivities on the performance of his ACAP technique, it was shown that even a considerable level of noise could be tolerated. The technique does not take into account the noise in the function values, i.e. deals with the noise in the design sensitivities only.

The performance of global approximation techniques was investigated by Schoofs (1987), Giunta et al. (1994), Narducci et al. (1995), Roux et al. (1996) and others, see reviews by Barthelemy and Haftka (1993) and Sobieszanski-Sobieski and Haftka (1997). These techniques are based on the least-squares surface fitting approach which was originally developed for problems where the function values are obtained as a result of physical measurements containing some level of noise. This approach allows to construct explicit approximations valid in the entire design variable space, but it is becoming computationally expensive as the number of design variables grows.

The issue of domain-dependent calculability arises when at some points of the design variable space the response functions cannot be evaluated by the adopted simulation code, Toropov et al. (1999). A traditional measure to prevent the failure of the optimization process when such a problem is encountered is to try to reduce the step length in the optimization technique. However, this may considerably slow down the optimization process and, still, does not guarantee the trouble-free convergence. Physically, the non-existence of the response function means that an unrealistic design has been considered which cannot be analysed by a response analysis procedure. Often, such situations lead to a premature termination or failure of the optimization method. Examples of such problems include optimization of dynamic behaviour of nonlinear mechanical systems, Markine (1999), optimization of aerodynamic characteristics of aerofoils, Toropov et al. (1999), to name but a few.

This study describes the use of the multipoint approximation method based on the response surface fitting (Toropov 1989, Toropov et al. 1993, 1995, 1996, van Keulen and Toropov 1997), sometimes abbreviated to MAM or MARS.

This technique replaces the original optimization problem by a succession of simpler mathematical programming problems. The functions in each iteration present mid-range approximations to the corresponding original functions. These functions are noise-free or, at least, the level of noise does not influence the convergence of an individual optimization sub-problem. The solution of an individual sub-problem becomes the starting point for the next step, the move limits are changed and the optimization is repeated iteratively until the optimum is reached. Each approximation function is defined as a function of design variables as well as a number of tuning parameters. The latter are determined by the weighted least squares surface fitting using the original function values (and their derivatives, when available) at several points of the design variable space. This selection of points will be referred to as a plan of numerical experiments.

In order to deal with the issue of domain-dependent calculability, the points of a plan where a response function is impossible to compute (e.g. resulting in a crash of a simulation code) are to be identified in the process of planning of numerical experiments, and a different point allocated.

There are several important aspects in the MARS which considerably influence the rate of convergence, namely:

Before the discussion of these aspects, the outline of the method is presented in Section 2.

The MARS by its very concept is well suited for the implementation in a parallel computer environment. The main advantage is that, even with a simple and straightforward implementation, good scalability of the optimization process can be achieved. This aspect will be the subject of Section 7.

References to various applications of MARS will be presented in Section 8, and a final discussion will be given in Section 9.


For self-containment the MARS is summarized here. More details on various aspects of the method will be given in the subsequent sections.

A general optimization problem can be formulated as


where x refers to the design variables. In order to reduce the number of calls for the response function evaluations and to lessen the influence of noise, the MARS replaces the optimization problem by a sequence of approximate optimization problems:


where k is the iteration number.

The selection of the noise-free approximate response functions is such that their evaluation is inexpensive as compared to the evaluation of the response functions Fj, although they are not necessarily explicit functions of the design variables. The approximate response functions are intended to be adequate in a current search sub-domain. This is achieved by appropriate planning of numerical experiments and use of move limits and . Moreover, it is attempted to achieve the best quality of the approximation functions in those regions of the design variable space where the solution of the approximate optimization problem can be expected, e.g. on the boundary of the feasible region.

The approximations are determined by means of weighted least squares, which reads


 Here minimization is carried out with respect to the tuning parameters aj and wpj refers to the weight coefficients. Their selection will be discussed in Section 4.

If the first order derivatives at plant points xp are known, the least-squares surface fitting is equivalent to the minimization of the function


where is a normalizing coefficient and 0 < g <1 reflects the weight of the information on sensitivities as compared to the information on the function values.



The approximate response functions must be flexible enough to mimic the behaviour of the response functions with sufficient accuracy. On the other hand, they must be easy to evaluate and possess only a minor level of numerical noise, if any.  

Intrinsically Linear Approximations

A simple but usually efficient choice of the structure of the approximations is an intrinsically linear (with respect to the tuning parameters) model, e.g. a linear and a multiplicative model:

 and (5)

 Intrinsically linear functions have been successfully used for a variety of design optimization problems by Toropov (1989), Toropov et al. (1993, 1996), Toropov and Carlsen (1994), Polynkine et al. (1995, 1996a, 1996b), Markine et al. (1996a, 1996b), van Keulen et al. (1997), van Keulen and Toropov (1997), Markine (1999). The advantage of these approximation functions is that a relatively small number of tuning parameters ai is to be determined and the corresponding least squares problem is solved easily. After application of an appropriate transformation only a set of linear equations has to be solved. It should be mentioned, though, that the transformation changes the definition of the weight coefficients wpj used in the weighted least squares fitting.

The use of linear (in x) approximations often leads to the deterioration of the convergence at the end of an optimization process. Generally, the multiplicative approximations show better convergence characteristics but the function of the optimization problem should be normalized, for otherwise a floating point overflows can be encountered. Higher order approximations, e.g. full quadratic polynomials, lead to a significantly better quality of approximation. However, these approximation functions require a considerably larger number of designs to be evaluated.

It should be noted that the requirements to the accuracy of the mid-range approximations can be more relaxed as compared to the case of the global ones but, obviously, more accurate approximations would produce a faster convergence. Several new approaches are being investigated in the attempt to produce new high quality approximations valid for a larger range of design variables.

 Use of Genetic Programming Methodology

Selection of the structure of an analytical approximation function is a problem of empirical model building (Box and Draper, 1987). Selection of individual regression components in a model results in solving a combinatorial optimization problem. Even if the bank of all regressors is established (which is a difficult problem on its own), the search through all possible combinations would result in prohibitive computational effort. Toropov and Alvarez (1998) attempted to develop and use a genetic programming (GP) methodology for the creation of an approximation function structure of the best possible quality, and use it within a multipoint (or global) approximation technique.

GP is a branch of genetic algorithms (GA). While a GA uses a string of numbers to represent the solution, the GP creates a population of computer programs with a tree structure. In our case of design optimization, a program represents an empirical model to be used for approximation of a response function.

The programs are composed of elements from a terminal set, which includes the design variables x, and a functional set of mathematical operators that generate the regression model: { +, *, /, x y, etc. }. The tuning parameters a are allocated to the program according to the basic algebraic rules.

The evolution of the programs is performed through the action of the genetic operators (reproduction, crossover, mutation and elite transfer) and the evaluation of the fitness function which is a measure of the approximation quality, see the web site of Luis Alvarez for more details.

 Mechanistic Models

The idea to endow the approximation function with some properties of the original implicit function to improve the quality of the approximation stems from the empirical model-building theory. Box and Draper (1987) showed that a so-called mechanistic model, i.e. the one that is built upon some knowledge about the system under investigation, provides better approximations than purely empirical ones. An example of application of such a model to a problem of material parameter identification (formulated as an optimization problem) can be found in (Toropov and van der Giessen 1993) where explicit non-linear approximations of response quantities for a solid bar specimen descend from a functional dependence on material parameters for a simpler tubular specimen. A difficulty in using such models is that they depend on the specific features of the problem and the researcher's experience.

It should be stated that the approximation need not necessarily be an explicit function. It could be an implicit one if some numerical procedure is involved in its formulation. The basic requirements to such a model can be summarized as:

 Simplified Numerical Models

A more general way of constructing high quality approximations is to obtain a simplified numerical model by simplifying the analysis model (e.g. by using a coarser finite element mesh discretization, a reduced number of the natural modes of the model in dynamic analysis, simpler geometry, boundary conditions, etc.). Such a model should inherit the most prominent features of the original one and the response analysis with this model should be computationally much less expensive than with the original model. In such a case a simplified numerical model provides a good basis for development of high quality approximations (Toropov and van der Giessen 1993, Toropov and Markine 1996, Toropov et al. 1997).

A simplified numerical model can be then used to build the approximation model:  


 where f(x) is the function presenting the structural response using the simplified model. The approximation based on the simplified model satisfies all the above mentioned requirements. It reflects the main properties of the original complex model, it is computationally less expensive and relatively noiseless. Depending on the way of introduction and the number of the tuning parameters in the simplified expression (6), the following three types of the approximations have been proposed (Toropov and Markine, 1996):

 1. As a simplest case, the approximation function is a linear or multiplicative function of two tuning parameters:


 The vector of the tuning parameters then consists of two elements: a = [a0, a1]T.

 2. Alternatively, the tuning parameters can be introduced in an explicit correction function C(x,a), which also depends on the design variables and can be developed in exactly the same way as analytical approximation functions described above. For example, using the linear and multiplicative models as the correction function, the following approximations can be built:

 or .

 The same approach has been used by Venkataraman et al. (1998) who referred to C(x,a) as correction response surface approximations. 

3. The simplified model, providing the basis for the approximation, depends in fact not only on the design variables x but also on the other parameters of the finite

element model such as material and/or geometry properties. To fit the simplified model into the data obtained with the complex model at the plan points, some of the parameters can be considered as the tuning parameters. The choice of the tuning parameters a should then be based on physical grounds. The approximation model can then be written in the form  




The weight coefficients influence the relative contribution of information on function values (and their derivatives, if available). Their choice strongly influences the quality of the approximations in different regions of the design variable space. Since the optimum point usually belongs to the boundary of the feasible region, the approximation functions should be more accurate in that region. Thus, the information at the points located near the boundary of the feasible region is to be treated with greater weights. This can be achieved by allocating an appropriate weight to the value of a constraint function at a point, e.g. .

The weights will then depend on the values of the constraint functions. It should be noted that in previous work some other definitions of weights have been applied., see, for example, Toropov et al. (1993).

In a similar manner a larger weight can be allocated to a design with a better objective function, see Toropov et al. (1995), van Keulen et al. (1996, 1997), van Keulen and Toropov (1997) for more details. Numerical examples showed, however, that the effect of this additional weight coefficient is relatively small.

In certain cases, e.g. when adaptive mesh refinement techniques are applied, it is possible to specify the error requested from the numerical analysis procedure (error of the numerical experiment). The difference in the quality of data, produced by response analyses with different prescribed errors, can also be taken into account by means of an additional weight coefficient . This can be efficiently used for the solution of problems where the response analysis is obtained using models of variable complexity and accuracy. Van Keulen et al. (1997) and van Keulen and Toropov (1997) showed examples where a good strategy for the choice of the prescribed error (i.e. complexity of the model) throughout the optimization process could significantly reduce the overall computational effort. The basic idea is not to keep a constant prescribed error (and model complexity) but to use a relatively rough model in the beginning of the optimization process and to refine it as the process progresses. Correspondingly, the quality of approximations can be low in the beginning and improve in the process of optimization.

The weight coefficients finally used for the weighted least squares fitting are now taken as .

When no information on the accuracy of the response functions is available or a certain contribution is not relevant, the corresponding term is taken as unity.

After a number of optimization steps have been carried out, a database with response function values becomes available. In addition to these design points available from the database, new designs are evaluated as part of a step in the optimization process. In order to achieve good quality approximations in a current search sub-domain, a selection of plan points must be made. All designs located in the current search sub-domain will be included in the weighted least squares fitting (if the corresponding response function has been successfully evaluated, of course). Points which are located far from the current search sub-domain would not generally contribute to the improvement of the quality of the resulting approximation functions in the current search sub-domain. For this reason only points located in the neighbourhood of the current search sub-domain are taken into account, as depicted in Fig. 7. A box in the space of design variables, which is approximately 1.5 to 2.0 times larger than the box representing the current search domain, was found by numerical experimentation to be a reasonable choice for the size of the neighbourhood.


Figure 7. Current search sub-domain and its neighbourhood (points marked ° are not included)

It should be noted that rejection of points not belonging to the neighbourhood could also be looked upon as if a zero weight was assigned to a particular design.



In the course of the optimization process new points must be generated in specified domains of the design variable space. In the literature several schemes have been proposed for generating new plans of experiments, see Burgee et al. (1996) among others. In previous work related to the MAM (Toropov 1989, Toropov et al. 1993 and other), a simple plan of experiments has been used which was based on n new points, where n equals the number of design variables. The locations of these points were obtained by a perturbation of each design variable by a fraction of the corresponding size of the current search sub-domain, as shown in Fig. 7 by larger circles. Generally, this scheme works well in conjunction with intrinsically linear approximations (5) but, due to its simplicity, it may have some disadvantages. Firstly, it does not take into account design points which are already located in the present search sub-domain and, therefore, newly evaluated designs may be located too close to previously evaluated designs. Secondly, the scheme is inflexible with respect to the number of added points and, therefore, has to be modified if other than intrinsically linear approximation functions are used.

In the present setting we intend to formulate a new scheme which does not suffer from the disadvantages mentioned above. The scheme must work efficiently even if only a single new design point is generated. Hence if there is a need for a certain number of new design points, the algorithm will be used several times, each time funding the same number of unknowns, namely n. Another reason for the development of a scheme adding a single point per run comes from the adopted parallel implementation, as will be discussed in Section 7. Finally, the scheme shall remain effective if different approximation functions are used within the same optimization task.

The simplest possible scheme is to add new points randomly. This is one of the ways of dealing with the issue of domain-dependent calculability of response functions (Toropov et al. 1999). Each added point is checked for calculability of the response function and, if the check fails, a new point is generated until a required number of plan points (all passing the check) are obtained. However, in this simple scheme the previously created points are not taken into account which can lead to poor plans.

The approach adopted here is to distribute points as homogeneously as possible in the sub-domains of interest. This is done by the introduction of a cost function

which is minimized with respect to the location of the new point d. Symbols denoted refer to coordinates which are normalized in the sub-domain of interest. The first term in the expression attempts to maximize the distance between points, and the second term promotes a homogeneous distribution along the coordinate axes. The third and fourth terms ensure that points do not belong to the boundary of the sub-domain. The last term prevents points from aligning along the diagonal of the search sub-region when only a few points are available. In the present implementation Q is minimized using a random search algorithm. We remind that if more points are to be allocated to a specific sub-domain, the minimization of Q is repeated to obtain coordinates for each new point.

As mentioned, in each step of the optimization a new set of designs, located in the current search sub-domain, will be evaluated. The number of additional points can be varied, e.g. in the current implementation the minimum number of new plan points equals the difference between the largest number of tuning parameters used in any of the approximation functions and the mumber of plan points already available in the current search sub-domain. As all points in the current search sub-domain and its neighborhood are involved in the weighted least squares fitting, more points are included in the plan of experiments than the minimum possible number. This has the advantage that the approximations become more reliable and it also reduces the effect of noise. Consequently, the convergence of the optimization process is improved.

Initial results indicate that it may be advantageous to add even more points during a step of the optimization process. These additional points can be located in the current search sub-domain, but may also be placed somewhere in the entire search domain. In the latter case the algorithm should be accompanied by a scanning process of the database. After solving the approximate optimization problem and evaluation of the corresponding design, the database is scanned for better points in terms of both objective and constraint functions. If a better point is found, this will serve as the starting point of the subsequent optimization step. If this happens, a jump in the optimization trajectory can be observed. The resulting algorithm thus could be looked upon as a blend of a random search and a successive approximation technique. Disadvantages are peculiar convergence characteristics and a higher cost of computations. The main advantage is that the algorithm tends to be more robust in terms of lesser sensitivity to local optima. So far, the question remains unanswered whether it is more efficient to add points outside the current search sub-domain as compared to a comprehensive initial scanning of the full search domain in order to find a good starting point; this is being currently investigated.

It should be stated that when design sensitivities are available, i.e. formulation (4) is used in the least-squares surface fitting, a plan of new numerical experiments becomes trivial: it adds the point obtained as the solution of the approximated problem (2) in a previous iteration to the data base of previously examined points.



After having solved the approximate optimization problem, a new search subregion must be defined, i.e. its new size and its location have to be specified. This is done on the basis of a move limit strategy which is summarized here. The reader is referred for more details to van Keulen and Toropov (1997) and van Keulen et al. (1996).  

Approximations "bad"



Stop: No convergence found



Reduce moderately or even enlarge


Approximation "reasonable"







Stop: Convergence found









Keep size

Approximations "good"







Stop: Convergence found








Keep Model & Enlarge

Keep Model & Keep size

 Table 1. Overview of move limit strategy

 The first step is to evaluate several indicators which are the basis for the move limit strategy. The first indicator is the quality index of the approximation functions. It is defined as the largest relative approximation error at the set of design variables corresponding to the solution of the approximate optimization problem. The quality of the approximation functions is then categorized as "bad", "reasonable" or "good". The second indicator is the location of the sub-optimum point in the current search subregion. When none of the current move limits is active, the solution is considered as "internal", otherwise it is denoted as "external". A third and fourth indicator are based on the movement history. For that purpose the angle between the last two move vectors is calculated. This will be the basis to identify the movement as "backward" or "forward". If the move vectors are nearly parallel, the convergence history will be labelled as "straight", otherwise it is marked as "curved". The fifth indicator, used in the termination criterion, is the size of the present search subregion. According to this indicator, the size can be "small" or "large". The sixth indicator is based on the value of the most active constraint. It is used to label the current solution as "close" or "far" from the boundary of the feasible and infeasible regions in the space of design variables.

Once the indicators have been determined, the sub-domain is moved and resized. In the present implementation only an isotropic resizing is used. Both reduction and enlargement of the sub-domain can be applied. Several levels of resizing can be applied, ranging from a moderate resizing, typically changing the size in the order of 10%, up to a quite aggressive resizing when the changes can be of the order of 100% or more. A summary of the move limit strategy as well as termination criteria is presented in Table 1.

If the obtained point does not pass the check for calculability of the response functions, the search sub-domain is reduced and the approximated problem (2) is solved again. 


In this a specific implementation of the MAM in a heterogeneous parallel computing environment is discussed.

As outlined in the previous sections, the algorithm includes two stages at which design evaluations are carried out: firstly, when new design points have been generated and, secondly, when the approximate nonlinear programming problem has been solved and a sub-optimum point obtained. At the former stage a relatively large number of designs are to be evaluated, whereas at the latter stage only a single design is evaluated. This single design evaluation is required to determine whether the iterative process is to be terminated and, if not, provides information required by the move limit strategy. The tasks of the first stage can be easily carried out in parallel, although a careful planning may be required to achieve a proper load balancing. The second phase is more difficult to execute in parallel as it requires a single design evaluation.

Fig. 8 illustrates a straightforward implementation without any planning or scheduling. It is seen that all available slaves (which can be processors of a multi-processor computer or computers in a network) are used to evaluate the new plan points. Thereafter the approximations are determined and the approximate optimization problem is solved. For the obtained design the response functions are evaluated by one of the slaves (the isolated block in the middle of Fig. 8). All remaining slaves are inactive at that time. It is seen that this simple implementation may lead to a relatively large idle time for individual nodes.


Figure 8. Example of a load on slaves in a heterogeneous network with five nodes. A block corresponds to the evaluation of a single design. Blue indicates an active slave, redindicates non-active slaves

 Apart from load balancing, the most obvious improvement would be the elimination of the single design evaluation which is carried out after identification and solution of the approximate optimization problem. However, this would make the move limit strategy much less effective as no assessment of the quality of approximations is then available. Instead, a different strategy for parallelization of the algorithm is suggested here. The idea is to use idle nodes for the enhancement of the database by adding selected response function evaluations as explained below. At the design evaluation stages of the optimization process an additional design evaluation is submitted to an idle node when at least two nodes become idle. The remaining idle node shall be kept ready for the execution of the sequential step of the algorithm, i.e. the design evaluation at the obtained sub-optimum point. The selection of coordinates of the additional points is carried out on the basis of the strategy described in Section 5. The size of the sub-domain to which the new points are to be added is taken identical to that of the current search sub-domain. Before solving the approximate optimization problem, the center of the sub-domain is taken as the best design evaluated so far. As soon as the approximate optimization problem has been solved, its solution is taken as the center. This rule leads to a good probability of putting the additional designs either in the next search sub-domain or sufficiently close to it. Thus, if an additional design belongs to the next search sub-domain, a smaller number of regular new design points have to be added to the plan corresponding to the next step of the optimization process. If an additional design falls outside the next search sub-domain, it may still contribute to the quality of the approximation functions, as all points belonging to the neighborhood of a search sub-domain are involved in the weighted least-squares fitting process. It should be noted that the proposed algorithm becomes most efficient when it is combined with the scanning mentioned in Section 5.

Figure 9 illustrates the proposed algorithm in a heterogeneous environment. All yellow blocks correspond to additional designs evaluated as free nodes become available. Obviously, the proposed algorithm could be combined with an appropriate load balancing technique.


Figure 9. Load on several slaves in a heterogeneous network with five nodes according to the proposed algorithm. A block corresponds to the evaluation of a single design. Blue indicates an active slave, red indicates non-active slaves. A yellow block corresponds to an active slave which is carrying out an additional design evaluation.



The performance of the present formulation of the MAM has been studied on a wide range of problems with noisy response analysis output. These studies include:


The multipoint approximation method, based on a sequence of approximate optimization problems, provides the possibility of dealing with noisy response functions effectively. This is to be attributed to the weighted least squares fitting used for determining the tuning parameters in the approximation functions. As soon as a sufficient number of design points, i.e. more than required minimum number, are included in this fitting, the method tends to average out noisy disturbances.

Settings for several coefficients in the algorithm have been suggested by van Keulen and Toropov (1997) on the basis of numerical experimentation with several simple optimization tasks with artificially introduced noise, and shape optimization tasks combined with AMR. More recent applications of MARS to shape optimization and optimization of thermoforming processes have not indicated any need for change of the coefficients determined earlier.

In contrast to previous work, a different scheme for the selection of new points in a plan of experiments has been adopted. It provides a flexible basis for either use of a small number of points in a plan or for the enhancement of the plan by additional points. The former is especially attractive if the level of noise is sufficiently small and many local optima are not expected. In such cases the use of a limited number of points is especially advantageous during the last few iterations of the optimization process. The previous formulation evaluates the same number of designs in each iteration which becomes too many near the optimal point. By using a limited number of new points, this effect can be eliminated.

The other possibility, i.e. adding more points than the required minimum number, has several effects. Firstly, better quality of the approximation functions can be achieved, especially when the response functions are significantly spoiled by the noise. Secondly, when it is combined with a scanning of the database, the chance for ending up in a local optimum can be reduced. When a problem has many local optimal points, it may even be advantageous to add a limited number of additional points in the entire search domain in each optimization step. Clearly, each step of the optimization process would then take longer time, but it could result in better convergence and the quality of the final design.

A non-standard implementation for a parallel computing environment has been selected. In this implementation idle slaves are used to evaluate designs which are expected to be close to the optimal point of a current step. Accordingly, these additional design evaluations will with a good probability contribute to the improvement of approximations used in the current and the next step. Moreover, as they are placed in regions of the design variable space which may contain better design alternatives than the present search sub-domain, these additional design points may produce the best design available in the database. When this is combined with a scanning of the database, the following observations have been made. First, the convergence history becomes somewhat unusual as it may now contain small jumps in the space of design variables, but the overall performance becomes better. Secondly, robustness is improved in terms of increasing the chance for finding a better non-local design solution.