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.

Collaboration:

burgeon

TU Delft

1. INTRODUCTION

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

the objective and constraint functions are evaluated as a result of
expensive numerical computations, e.g. using a FEM or a CFD code;

function values and/or their derivatives may contain some level of
numerical noise;

domain-dependent calculability of some of the response functions, i.e.
situations when these functions cannot be evaluated at some points of the
design variable space.

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 et.al. (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:

the choice of the structure of the approximate expressions (Section
3);

the strategy for allocation of the weight coefficients (Section 4);

the choice of the plan (design) of numerical experiments (Section 5);

the move limit strategy (Section 6).

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.

2. MULTIPOINT APPROXIMATION METHOD

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

(1)

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:

(2)

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

(3)

Here minimization is carried out with respect to the tuning parameters
ajand 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

(4)

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.

3. APPROXIMATIONS

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:

it must depend on the same design variables as the original function;

it has to contain some tuning parameters to be defined using the general
(non-linear) least-squares method;

it must be simple enough to be used in numerous repeated calculations;

it should not contain any considerable level of numerical noise in
order not to cause convergence problems in the optimization process.

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:

(6)

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:

or

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

.

4. WEIGHT COEFFICIENTS

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.

5. PLAN OF EXPERIMENTS

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.

6. MOVE LIMIT STRATEGY

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"

"Small"

"Large"

Stop: No convergence found

"Straight"

"Curved"

Reduce moderately or even enlarge

Reduce

Approximation "reasonable"

"Internal"

"External"

"Small"

"Large"

"Backward"

"Forward"

Stop: Convergence found

"Close"

"Far"

Reduce

"Straight"

"Curved"

Reduce

Reduce

Enlarge

Keep size

Approximations "good"

"Internal"

"Boundary"

"Small"

"Large"

"Backward"

"Forward"

Stop: Convergence found

"Close"

"Far"

Reduce

"Straight"

"Curved"

Reduce

Reduce

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.

7. PARALLEL IMPLEMENTATION

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.

8. APPLICATIONS

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:

Shape
optimization with AMR by van Keulen et al. (1997), van Keulen and Toropov
(1997, 1998), Polynkine et al. (1996). In these studies the accuracy of
the response analyses is controlled by the optimization process and, therefore,
is gradually improved in the course of optimization. This approach required
significant modifications to the move limit strategy and the selection
of the weights.

Optimization of products made of a continuous fibre reinforced thermoplastic
material by Polynkine et al. (1996). Here a sequence of numerical
simulations in the response analysis leads to a very noisy behaviour of
the response functions. Due to the high level of noise previous formulations
were easily trapped by local solutions. The present formulation was shown
to provide much better convergence characteristics. On a matter not directly
related to the optimization technique, it should be noted that some recent
improvements to the simulation schemes have reduced the level of noise.

Studies with artificially added noise by van Keulen and Toropov (1997)
and Toropov et al. (1996). These studies used an analytical noise-free
problem in which random noise was artificially added to the response functions.
A large number of optimization tasks has been performed in order to provide
an objective test bed for the selection of coefficients used in the algorithm.

Optimization
of aerofoil by Toropov et al. (1999) where, apart from dealing
with numerical noise, it was necessary to address the issue of domain-dependent
calculability of the objective and constraint functions produced by a CFD
code.

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.