# Parallel Surrogate-assisted Optimization Using Mesh Adaptive Direct Search

Bastien Talgorn1, Stéphane Alarie2, and Michael Kokkolaras1

(  1McGill University, GERAD, Montréal, Québec, Canada
2Hydro-Québec’s Research Institute, GERAD, Montréal, Québec, Canada
)

###### Abstract

We consider computationally expensive blackbox optimization problems and
present a method that employs surrogate models and concurrent computing at the
search step of the mesh adaptive direct search (MADS) algorithm. Specifically,
we solve a surrogate optimization problem using locally weighted scatterplot
smoothing (LOWESS) models to find promising candidate points to be evaluated
by the blackboxes. We consider several methods for selecting promising points
from a large number of points. We conduct numerical experiments to assess the
performance of the modified MADS algorithm with respect to available CPU
resources by means of five engineering design problems.

##  1 Introduction

We consider the optimization problem

$\begin{array}[]{rl}\underset{\mathbf{x}\in\mathcal{X}}{\min}&f(\mathbf{x})\\\
\text{subject to}&c_{j}(\mathbf{x})\leq 0,~{}~{}j=1,2,\ldots,m,\end{array}$
($P$)

where $f(\mathbf{x})$ is the objective function,
$\mathbf{x}\in{\mathbb{R}}^{n}$ is the vector of decision variables,
$\mathcal{X}$ is a subset of ${\mathbb{R}}^{n}$, and $c_{j}(\mathbf{x})$ are
general nonlinear constraints. We assume that some (at least one) of the
functions $\\{f,c_{1},c_{2},\ldots,c_{m}\\}$ are evaluated using simulations
or other computational procedures that are blackboxes. In particular, we
consider the case where these blackboxes are computationally expensive,
possibly nonsmooth and/or nonconvex, and that the process used to evaluate
them may crash or fail to return a value. Finally, we assume that function
gradients either do not exist theoretically or, if they do, cannot be computed
or approximated with reasonable computational effort.

Metaheuristics and derivative-free search algorithms are commonly used for
solving ($P$). The former (e.g., genetic algorithms(GAs), particle swarm
optimization (PSO), tabu search (TS), etc.) are commonly used for global
exploration while the latter (e.g., generalized pattern search (GPS), mesh
adaptive direct search (MADS), and trust-region methods (DFO, COBYLA, CONDOR))
are local methods with convergence properties [1]. In this work, we use the
NOMAD implementation [2, 3] of the MADS algorithm [4] to solve ($P$).

Multiprocessor computers, supercomputers, cloud computing, or just a few
connected PCs can provide parallel (or concurrent) computing opportunities to
speed up so-called trajectory-based optimization algorithms. According to [5],
three ways are commonly used to achieve this: (i) parallel evaluation of
neighborhood solutions (distributed evaluations), (ii) parallel trajectories
from the same (or different) initial guess(es) (independent optimization
runs), (iii) the evaluation of a point $\mathbf{x}$ is performed in parallel
(i.e., the search is sequential). The implementation of (iii) depends only on
the blackbox, while the other two are related to the optimization algorithm.

NOMAD offers implementations for (i) and (ii) through p-MADS for (i) and Coop-
MADS for (ii) [6]. In both cases, the parallel evaluations are handled by
NOMAD by means of MPI calls to the blackbox [3]. However, if one prefers to
take charge of the distribution of evaluations, one can implement p-MADS with
blocks by using NOMAD in batch mode [6]. In that case, instead of using MPI
calls, NOMAD writes all points to be evaluated in an input file, waits for all
evaluations to be completed, and reads the obtained values of $f(\mathbf{x})$
and $c_{j}(\mathbf{x})$ from an output file. Then, NOMAD either stops if a
local optimum is reached or submits a new block of points to evaluate. We use
here the block functionality of NOMAD, adopting option (i) for parallel
evaluations, because of its flexibility and generality.

The MADS algorithm includes two steps at each iteration, the SEARCH and the
POLL. The SEARCH step is flexible (defined by the user) and aims at
determining one or more new points $\mathbf{x}\in\mathcal{X}$ that improves
the current best solution. The POLL step is defined according to the
convergence analysis of the algorithm and generates trial points around the
current best solution. The number of poll points at each iteration is either
$2n$ or $n+1$ depending on the utilized pattern with $n$ being the size of
$\mathbf{x}$.

The number of points evaluated in the SEARCH step depends on the methods
chosen or defined by the user. Several techniques are already implemented and
available in NOMAD, including the speculative search (SS), the quadratic model
search (QUAD), the variable neighborhood search (VNS), the Latin hypercube
search (LHS), and the Nelder Mead search (NMS). One can also implement their
own technique if one so desires, which is called a _user search_ (US).

With the exception of LHS, all provided techniques usually return only one
trial point. When several techniques are used at once, they are called one
after the other along the SEARCH step, each technique providing its own trial
point, which is evaluated by the blackbox before proceeding to the next
technique. Assuming that $2n$ CPUs are available for solving ($P$), the POLL
step can make good use of these CPUs. However, since SEARCH step evaluations
are sequential, progress is slowed down with almost all CPUs being idle. One
may argue that we should then only use LHS since it can generate $2n$ points.
However, since LHS is random, its points will quickly become less promising
after a few iterations.

Considering that the number of available CPUs are now, particularly with the
emergence of cloud computing, relatively inexpensive and unlimited, we should
rethink the SEARCH step to be as effective as the POLL step in terms of CPU
use. In this work, we propose a SEARCH step technique that returns a large
number of diverse points for evaluation.

The paper is structured as follows. The general idea behind the proposed
technique is described in Section 2. In Section 3, six different methods are
presented for selecting various candidates from a large set of points. In
Section 4, the resulting model search for parallel computing is specified. In
Section 5, we test our SEARCH step technique on five engineering design
optimization problems using up to 64 processors. A discussion concludes the
paper.

##  2 Proposed SEACH step technique

One of the practical challenges of the SEARCH step is that only one candidate
is obtained at a significant computational investment [7, 8, 9, 10].
Specifically, regardless of the number of available CPUs, only one CPU is used
in the SEARCH step for blackbox evaluations, with the exception of LHS. Before
presenting our idea for mitigating this practical challenge, we will assume
that computationally inexpensive surrogate models of the expensive blackboxes
are available. We can then consider the surrogate problem of problem ($P$)

$\begin{array}[]{rl}\underset{\mathbf{x}\in\mathcal{X}}{\min}&\hat{f}(\mathbf{x})\\\
\text{subject to}&\hat{c}_{j}(\mathbf{x})\leq
0,~{}~{}j=1,2,\ldots,m,\end{array}$ ($\hat{P}$)

where $\\{\hat{f},\hat{c}_{1},\hat{c}_{2},\ldots,\hat{c}_{m}\\}$ are surrogate
models of $\\{f,c_{1},c_{2},\ldots,c_{m}\\}$, respectively. Note that we only
need to ensure that the minimizers of ($P$) and ($\hat{P}$) are close enough,
and not that the surrogate models are good approximations of the blackboxes
globally. It then follows that a minimizer of ($\hat{P}$) will be a good
candidate for the solution of ($P$).

If both problems have the same minimizers, they may share features in other
areas of $\mathcal{X}$ as well. Since the evaluations of $\hat{f}(\mathbf{x})$
and $\hat{c}_{j}(\mathbf{x})$ are rather inexpensive compared to $f$ and
$c_{j}$, one can allow a very large budget of model evaluations to solve
($\hat{P}$), extending thus the number of design space areas that will be
visited. This is acceptable as long as the solution of ($\hat{P}$) is faster
than any single evaluation of the blackboxes. Considering there are $q$ CPUs
available for blackbox evaluations, one may then select $q$ points from the
available budget by solving ($\hat{P}$). The $q$ points can be selected to
consider areas of $\mathcal{X}$ that have been neglected until now in the
solution of ($P$).

The above proposition proposes the use of surrogate models
$\\{\hat{f},\hat{c}_{1},\hat{c}_{2},\ldots,\hat{c}_{m}\\}$ in a manner that is
not reported in [11], which mentions two ways of exploiting surrogates in the
context of parallelization. The simplest is to fit $q$ different surrogate
models at the same points already evaluated by the blackbox functions. This
allows to get $q$ different promising candidates and requires no uncertainty
quantification for the surrogate models. One can also combine the surrogates;
distance-based criteria can be added to ensure diversity between the
candidates. The other way is to use a single surrogate model and consider $q$
points where the blackboxes should be evaluated at to improve its accuracy.

We propose an intermediate approach. We use only one surrogate model for each
blackbox. The $q$ candidates are extracted from that single surrogate, but not
with the aim of improving it. Instead, the $q$ candidates are selected to be
the most interesting to advance the optimization process.

##  3 Methods for selecting candidate points

Let $\mathbf{X}=\\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{k}\\}$ be
the set of all points evaluated by the blackbox. Note that
$\mathbf{X}\subset\mathcal{X}\subset{\mathbb{R}}^{n}$. We denote
$\mathbf{\hat{X}}$ the _surrogate cache_ , i.e., the set of all points for
which $\\{\hat{f},\hat{c}_{1},\hat{c}_{2},\ldots,\hat{c}_{m}\\}$ have been
evaluated during the solution of ($\hat{P}$). Similarly,
$\mathbf{\hat{X}}\subset\mathcal{X}\subset{\mathbb{R}}^{n}$. Let $\mathbf{S}$
be the set of points that are selected by the SEARCH step to be evaluated with
the blackbox. The set $\mathbf{S}$ is initially empty and is built from the
points of $\mathbf{\hat{X}}$ (ensuring that
$\mathbf{S}\subset\mathbf{\hat{X}}$) with a greedy algorithm by means of up to
six selection methods, each one having a different goal.

  * •

Method 1 selects the best point of $\mathbf{\hat{X}}$ not in
$\mathbf{X}\cup\mathbf{S}$;

  * •

Method 2 selects the most distant point of $\mathbf{\hat{X}}$ from
$\mathbf{X}\cup\mathbf{S}$;

  * •

Method 3 selects the best point of $\mathbf{\hat{X}}$ at a certain distance of
$\mathbf{X}\cup\mathbf{S}$;

  * •

Method 4 selects the best point of $\mathbf{\hat{X}}$ under additional
constraints;

  * •

Method 5 selects a point of $\mathbf{\hat{X}}$ that is a possible local
minimum of the surrogate problem;

  * •

Method 6 selects a point of $\mathbf{\hat{X}}$ in a non-explored area.

Note that some selection methods may fail to return a candidate, particularly
methods 3 and 4. If this happens, the next method is used. We repeat and loop
through all methods until we obtain enough points in $\mathbf{S}$ matching the
available CPUs. Some points $\mathbf{x}\in\mathbf{\hat{X}}$ can also belong to
the set $\mathbf{X}$; this is not an issue since all methods only select point
$\mathbf{x}\in\mathbf{\hat{X}}$ to be added to $\mathbf{S}$ if and only if
$\mathbf{x}\notin\mathbf{X}$. The selection methods are detailed below after
some definitions.

###  3.1 Definitions

Let $d(A,B)$ be the Euclidean distance between two subsets $A$ and $B$ of
${\mathbb{R}}^{n}$

$d(A,B)=\underset{a\in A}{\min}\;\underset{b\in B}{\min}\;\|a-b\|_{2}.$ (1)

As a convention, the distance to an empty set is infinite:
$d(A,\varnothing)=d(\varnothing,\varnothing)=+\infty$. By extension, we will
denote the distance between an element $a\notin B$ and the subset $B$ simply
by $d(a,B)$, which implies that $a$ also refers to the particular subset
containing only $a$, i.e., $\\{a\\}$.

Regarding feasibility, we consider the aggregate constraint violation function
used in [12], i.e.,
$h(\mathbf{x})=\sum_{j=1}^{m}\max\\{0,c_{j}(\mathbf{x})\\}^{2}$. The same
function is used in the _progressive barrier_ mechanism in NOMAD [13].

We also define the order operators between two points $\mathbf{x}$ and
$\mathbf{x}^{\prime}\in\mathcal{X}$:

$\displaystyle\mathbf{x}\prec\mathbf{x}^{\prime}\Leftrightarrow$
$\displaystyle\left\\{\begin{array}[]{l}h(\mathbf{x})<h(\mathbf{x}^{\prime})\\\
\text{or}\\\ h(\mathbf{x})=h(\mathbf{x}^{\prime})\text{ and
}f(\mathbf{x})<f(\mathbf{x}^{\prime})\end{array}\right.$ (5)
$\displaystyle\mathbf{x}\preceq\mathbf{x}^{\prime}\Leftrightarrow$
$\displaystyle\;\;{\tt not}(\mathbf{x}^{\prime}\prec\mathbf{x}),$ (6)

which are transitive. By those definitions, an incumbent solution
$\mathbf{x}^{\prime}$ of the original problem ($P$) is such that
$\mathbf{x}^{\prime}\preceq\mathbf{x},\;\forall\mathbf{x}\in\mathbf{X}$.
Similarly, a global minimizer $\mathbf{x}^{*}$ is such that
$\mathbf{x}^{*}\preceq\mathbf{x},\;\forall\mathbf{x}\in\mathcal{X}$. In the
same manner as we define $\prec$ and $\preceq$ for $f$ and $h$, we define
$\;\widehat{\prec}\;$ and $\;\widehat{\preceq}\;$ for $\hat{f}$ and $\hat{h}$.

Finally, to simplify the description of the proposed selection methods, we
define $\mathbf{s}^{\infty}$ as a _virtual point_ (in the sense that it does
not have coordinates in ${\mathbb{R}}^{n}$), which represents the worst
possible candidate in ${\mathbb{R}}^{n}$:

$\hat{h}(\mathbf{s}^{\infty})=\hat{f}(\mathbf{s}^{\infty})=+\infty\text{ \;\;
and \;\; }d(\mathbf{s}^{\infty},\mathbf{X})=0.$ (7)

###  3.2 Detailed description of selection methods

Method 1. The first selection method selects the best point $\mathbf{s}$ of
$\mathbf{\hat{X}}$ under the constraint that
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})>0$, which means that $\mathbf{s}$ is
not in the set $\mathbf{X}$ of evaluated points nor already selected (i.e.,
$\notin\mathbf{S}$). This method reflects how surrogate models are typically
used for finding new candidate points.

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{for all $\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\
\left|\;\begin{array}[]{l}\text{if
$\mathbf{s}\;\widehat{\prec}\;\mathbf{s}^{*}$ ~{}{\tt and}~{}
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})>0$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}$

Algorithm 1 Selection of the best point (Method 1)

Method 2. The second method aims to maximize the diversity of the candidates
to be evaluated. It selects the point $\mathbf{s}$ of $\mathbf{\hat{X}}$ that
maximizes the distance $d(\mathbf{s},\mathbf{X}\cup\mathbf{S})$, i.e., as far
as possible from points already evaluated.

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{for all $\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\
\left|\;\begin{array}[]{l}\text{if
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})>d(\mathbf{s}^{*},\mathbf{X}\cup\mathbf{S})$,
then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}$

Algorithm 2 Selection of the most distant point to $\mathbf{X}\cup\mathbf{S}$
(Method 2)

Method 3. This method selects the best point $\mathbf{s}$ of
$\mathbf{\hat{X}}$ under the constraint that
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})\geq d_{\min}$, where $d_{\min}$ is
initialized at 0 at the beginning of the selection process and increased
progressively as the method is applied. Method 3 may fail to select a
candidate $\mathbf{s}$ when $d_{\min}$ becomes too large. Since the selected
points $\mathbf{S}$ must be projected on the current mesh
$\mathcal{M}=\\{\mathbf{x}+\Delta^{m}\mathbf{D}\mathbf{z},\mathbf{z}\in{\mathbb{N}}^{n_{D}},\mathbf{x}\in\mathbf{X}\\}$
as required by MADS, incrementing $d_{\min}$ by the current mesh size
$\Delta^{\mathcal{M}}$ allows to avoid that several candidates become
identical after the projection.

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{if first use of Method~{}\ref{algo:method3}, then:}\\\
\left|\;\begin{array}[]{l}\text{$d_{\min}\leftarrow 0$}\\\
\end{array}\right.\\\ \text{end}\\\ \text{for all
$\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\ \left|\;\begin{array}[]{l}\text{if
$\mathbf{s}\;\widehat{\prec}\;\mathbf{s}^{*}$ ~{}{\tt and}~{}
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})\geq d_{\min}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\text{$d_{\min}\leftarrow d_{\min}+\Delta^{\mathcal{M}}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}$

Algorithm 3 Selection of the best point with a constraint on the distance to
$\mathbf{X}\cup\mathbf{S}$ (Method 3)

Method 4. Considering that the surrogate models $\hat{c}_{j}$ may fail to
predict correctly if $\mathbf{s}$ is feasible, the present method tries to
select points that will be likely to be feasible when evaluated by the
blackboxes $c_{j}$. This is done by selecting the best feasible point
$\mathbf{s}$ of $\mathbf{\hat{X}}$ under the constraint
$\hat{c}_{\max}(\mathbf{s})\leq\hat{c}_{\text{margin}}$, where
$\hat{c}_{\max}(\mathbf{s})$ is defined as being the most violated constraint
of $\mathbf{s}$, i.e.,

$\hat{c}_{\max}(\mathbf{s})=\underset{j=1,2,\dots,m}{\max}\;\hat{c}_{j}(\mathbf{s}),$
(8)

where $\hat{c}_{\text{margin}}$ is set as

$\hat{c}_{\text{margin}}\leftarrow\underset{\begin{subarray}{c}\mathbf{s}\in\mathbf{\hat{X}}\\\
\hat{c}_{\max}(\mathbf{s})<0\end{subarray}}{\max}\hat{c}_{\max}(\mathbf{s})$
(9)

and quantifies, among all feasible points of $\mathbf{\hat{X}}$, the smallest
amount by which these are satisfied.

By definition, $\hat{c}_{\max}(\mathbf{s})\leq 0$ if $\mathbf{s}$ is predicted
to be feasible by the surrogate models. The more negative
$\hat{c}_{\max}(\mathbf{s})$ is, the more likely is $\mathbf{s}$ to be
feasible when evaluated by the blackboxes. Decreasing progressively the value
of $\hat{c}_{\text{margin}}$ after each call of this selection method will
favor candidates that are increasingly likely to be feasible (but possibly
with worse objective function values).

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{if first use of Method~{}\ref{algo:method4}, then:}\\\
\left|\;\begin{array}[]{l}\text{$\hat{c}_{\text{margin}}\leftarrow\min\\{0,\underset{\begin{subarray}{c}\mathbf{s}\in\mathbf{\hat{X}}\\\
\hat{c}_{\max}(\mathbf{s})<0\end{subarray}}{\max}\hat{c}_{\max}(\mathbf{s})\\}$}\\\
\end{array}\right.\\\ \text{end}\\\ \text{for all
$\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\ \left|\;\begin{array}[]{l}\text{if
$\hat{c}_{\max}(\mathbf{s})\leq\hat{c}_{\text{margin}}$ ~{}{\tt and}~{}
$\hat{f}(\mathbf{s})<\hat{f}(\mathbf{s}^{*})$ ~{}{\tt and}~{}
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})>\Delta^{\mathcal{M}}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\text{$\hat{c}_{\text{margin}}\leftarrow
2\,\hat{c}_{\max}(\mathbf{s}^{*})$}\\\ \end{array}\right.\\\ \text{end}\\\
\end{array}$

Algorithm 4 Selection of the best point with a constraint on the feasibility
(Method 4)

Note that this method requires that $\hat{c}_{\text{margin}}\leq 0$ is always
satisfied. Moreover, we also assume that there is at least one
$\mathbf{s}\in\mathbf{\hat{X}}$ that is feasible. If it is not the case, which
may happen in the first iteration, we will end up with
$\hat{c}_{\text{margin}}>0$ and an inappropriate candidate will be selected.
To avoid this, we initialize $\hat{c}_{\text{margin}}$ to $0$ so that the
method may fail to return a candidate if that is the case.

Method 5. The isolation distance is used here to detect local minima of the
surrogate problem. This concept is inspired from the _topographic isolation_
of a mountain summit, which measures the local significance of a summit. It is
defined as the distance to the closest higher summit.111 In mountaineering,
the topographic isolation of Mount Everest is infinite and the summit with the
second highest isolation is the Aconcagua in Argentina. The Aconcagua is not
the second highest summit but there is no higher mountain in a 16,518 km
range, making it the most important summit in the Americas and in the southern
hemisphere.

Transferred to optimization, the concept of topographic isolation is used to
quantify the importance of a local minimum. Its strict application is however
impossible since it will require to prove that no other point within a certain
distance of $\mathbf{x}$ is better than $\mathbf{x}$. We can only compute
isolation distance of the already evaluated points. Consequently, we define
the isolation distance as being the distance from $\mathbf{s}$ to the closest
point of $\mathbf{\hat{X}}$ that is better than $\mathbf{s}$

$d_{\text{iso}}(\mathbf{s})=\underset{\begin{subarray}{c}\mathbf{s}^{\prime}\in\mathbf{\hat{X}}\\\
\mathbf{s}^{\prime}\\!\;\widehat{\prec}\;\\!\mathbf{s}\end{subarray}}{\min}\;d\big{(}\mathbf{s},\mathbf{s}^{\prime}).$
(10)

Constraints are taken into account by using the order relationship defined in
Equation (5). As a convention, if no point of $\mathbf{\hat{X}}$ is better
than $\mathbf{s}$, then $d_{\text{iso}}(\mathbf{s})=+\infty$. With this
definition, the point of $\mathbf{\hat{X}}$ with the highest isolation
distance is also the best candidate in $\mathbf{\hat{X}}$. However, we have
observed that the other points with a high isolation distance are often poor
points far from any other point of $\mathbf{\hat{X}}$. To address this
problem, we define the _isolation number_ of $\mathbf{s}\in\mathbf{\hat{X}}$
as the number of points of $\mathbf{\hat{X}}$ within the ball of centre
$\mathbf{s}$ and radius $d_{\text{iso}}(\mathbf{s})$

$n_{\text{iso}}(\mathbf{s})=\text{card}\big{\\{}\mathbf{s}^{\prime}:\mathbf{s}^{\prime}\in\mathbf{\hat{X}},d(\mathbf{s},\mathbf{s}^{\prime})<d_{\text{iso}}(\mathbf{s})\big{\\}}.$
(11)

To have a high isolation number, a point must be better than many of its
neighbors, which means that this criterion allows to detect local minima. Note
that Equation (7) implies that
$d_{\text{iso}}(\mathbf{s}^{\infty})=n_{\text{iso}}(\mathbf{s}^{\infty})=0$.
Method 5 selects the point of $\mathbf{\hat{X}}$ that has the highest
isolation number not yet in $\mathbf{X}\cup\mathbf{S}$.

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{for all $\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\
\left|\;\begin{array}[]{l}\text{if
$n_{\text{iso}}(\mathbf{s})>n_{\text{iso}}(\mathbf{s}^{*})$ ~{}{\tt and}~{}
$d(\mathbf{s},\mathbf{X}\cup\mathbf{S})>0$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}$

Algorithm 5 Selection of the most isolated point (Method 5)

Method 6. The purpose of this method is to select points in neglected areas of
the design space. To do so, it selects points in areas heavily explored while
solving ($\hat{P}$) but overlooked when solving ($P$). The _density number_ of
$\mathbf{s}\in\mathbf{\hat{X}}$ is defined as

$n_{\text{density}}(\mathbf{s})=\text{card}\big{\\{}\mathbf{s}^{\prime}:\mathbf{s}^{\prime}\in\mathbf{\hat{X}},d(\mathbf{s},\mathbf{s}^{\prime})<d(\mathbf{s},\mathbf{X}\cup\mathbf{S})\big{\\}}.$
(12)

Method 6 selects the point of $\mathbf{\hat{X}}$ with the highest density
number. Note that, as for $n_{\text{iso}}$, Equation (7) implies that
$n_{\text{density}}(\mathbf{s}^{\infty})=0$.

$\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}^{\infty}$}\\\
\text{for all $\mathbf{s}\in\mathbf{\hat{X}}$, do:}\\\
\left|\;\begin{array}[]{l}\text{if
$n_{\text{density}}(\mathbf{s})>n_{\text{density}}(\mathbf{s}^{*})$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{s}^{*}\leftarrow\mathbf{s}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}\right.\\\ \text{end}\\\
\text{if $\mathbf{s}^{*}\neq\mathbf{s}^{\infty}$, then:}\\\
\left|\;\begin{array}[]{l}\text{$\mathbf{S}\leftarrow\mathbf{S}\cup\\{\mathbf{s}^{*}\\}$}\\\
\end{array}\right.\\\ \text{end}\\\ \end{array}$

Algorithm 6 Selection of a point in a populated area (Method 6)

##  4 Parallel computing implementation

We now describe how we implement the proposed SEARCH step. We start with the
surrogate models and follow with the algorithms used for solving ($\hat{P}$).
We conclude with how all of this is integrated with the MADS algorithm to
solve ($P$).

###  4.1 Surrogate models

Several surrogate models are mentioned in [11] regarding their use in parallel
computing approaches, including Kriging, radial basis functions (RBF), support
vector regression (SVR), and polynomial response surfaces (PRSs). For a more
exhaustive description and comparison of surrogate models, see [14]. Based on
our previous work reported in [10], we choose to use the locally weighted
scatterplot smoothing (LOWESS) surrogate modeling approach [15, 16, 17, 18].

LOWESS models generalize PRSs and kernel smoothing (KS) models. PRSs are good
for small problems, but their efficacy decreases for larger, highly nonlinear,
or discrete problems. KS models tend to overestimate low function values and
underestimate high ones, but usually predict correctly which of two points
yields the smallest function value [9]. LOWESS models build a linear
regression of kernel functions around the point to estimate. They have been
shown to be suitable for surrogate-based optimization [10]. Their parameters
are chosen using an error metric called “aggregate order error with cross-
validation” (AOECV), which favors equivalence between the original problem
($P$) and the surrogate problem ($\hat{P}$) [9]. We use the SGTELIB
implementation of the surrogate models, which is now integrated as a surrogate
library in NOMAD version 3.8 [19]. Specifically, the considered LOWESS model
is defined in SGTELIB as follows.

$\text{\tt TYPE Lowess DEGREE 1 RIDGE 0 SHAPE\\_COEF OPTIM KERNEL\\_TYPE
OPTIM},$

which means that the local regression is linear, the ridge (regularization)
coefficient is 0, and the kernel shape and kernel type are optimized to
minimize the aggregate order error (see [10]).

The LOWESS model is built as described in Appendix A. Only the Gaussian kernel
was considered in [10]. Six additional kernel functions have meanwhile been
implemented in SGTELIB. Accordingly, not only $\lambda$ (kernel shape) is
chosen to minimize AOECV, but also $\phi$ (kernel type).

###  4.2 Surrogate problem solution

The surrogate problem ($\hat{P}$) is solved by means of an inner instance of
MADS; it is initialized by a Latin hypercube search (LHS) [20, 21] and uses
variable neighborhood search (VNS) [22, 23] as the SEARCH step and a large
budget of function evaluations (10,000). The POLL step is performed using the
ORTHO 2N directions option. This inner MADS is implemented in the SEARCH step
of the outer MADS.

The LHS guarantees that there are surrogate cache points widely spread over
the design space. To ensure this, 30% of all function evaluations (i.e.,
3,000) are devoted to the LHS. Four additional points are considered: The
current best feasible point of the original problem, the current best
infeasible point of the original problem, the best feasible point obtained by
solving the most recent surrogate problem instantiation, and the best
infeasible point by solving the most recent surrogate problem instantiation.
These points are used as initial guesses of the inner MADS problem, which will
be run until the remaining evaluation budget is exhausted. This budget will be
shared between the POLL step and the VNS in a default proportion where 75% is
devoted to VNS, which favors the exploration of multiple local attraction
basins. A large number of evaluations that build the surrogate cache
$\mathbf{\hat{X}}$ favors an accurate solution of the surrogate problem. Using
LHS, VNS, and a large number of function evaluations ensures that
$\mathbf{\hat{X}}$ contains highly promising candidates for the solution of
the original problem ($P$).

###  4.3 The modified MADS algorithm

Recall that each iteration of MADS includes a SEARCH step (performed first)
and a POLL step. Let $t$ denote the current iteration. Then, $\mathbf{X}_{t}$,
$\mathbf{\hat{X}}_{t}$, and $\mathbf{S}_{t}$ denote the sets $\mathbf{X}$,
$\mathbf{\hat{X}}$, and $\mathbf{S}$ considered by MADS at iteration $t$. Let
$q$ be the number of available CPUs, i.e., the number of blackbox evaluations
that can be performed in parallel. The proposed MADS for exploiting $q$ CPUs
is as follows. First, the SEARCH step proceeds by solving the surrogate
problem to populate the set $\mathbf{\hat{X}}_{t}$. From that set, $q$
candidates are selected and returned to be evaluated by the blackbox(es) in
parallel. The selection is made by cycling through a user-defined subset of
the six proposed selection methods (Section 3.2) until a total of $q$
candidates are selected, or until all selection methods consecutively failed
to add a candidate to $\mathbf{S}_{t}$. If $q$ is smaller than the number of
selection methods retained by the user, we do not necessarily go through all
the methods, but stop as soon as we get $q$ candidates.

If/when the SEARCH step fails to return a better objective function value, the
MADS algorithm proceeds to the POLL step. Let $\mathbf{P}_{t}$ be the set of
candidates produced by the polling directions at the iteration $t$. The
cardinality of $\mathbf{P}_{t}$ is denoted by $|\mathbf{P}_{t}|$. To be
consistent with our need to fulfill continuously all the available CPUs with
evaluations, additional candidates are added so that $|\mathbf{P}_{t}|$ is at
least $q$ or a multiple of $q$. This is accomplished by means of NOMAD’s
intensification mechanism ORTHO 1. If $|\mathbf{P}_{t}|=q$, all poll
candidates of $\mathbf{P}_{t}$ are evaluated concurrently, eliminating the
need to order them. If $|\mathbf{P}_{t}|>q$, then the points in
$\mathbf{P}_{t}$ are regrouped in several blocks of $q$ candidates. The blocks
are then evaluated sequentially and opportunistically, which means that if a
block evaluation leads to a success, the remaining blocks are not evaluated.
To increase the probability of success from the first block, and hence
avoiding to proceed with the remaining ones, the candidates of
$\mathbf{P}_{t}$ are sorted using the surrogate models and distributed in the
blocks so that the more promising ones are in the first block.

Recall that $\mathcal{M}_{t}$ is the current mesh, $\Delta^{\mathcal{M}}_{t}$
is the associated mesh size parameter, $\Delta^{\text{P}}_{t}$ is the
corresponding mesh poll parameter, and $\mathbf{x}^{*}_{t}$ is the best
solution found at iteration $t$. Finally, the set $\mathbf{X}_{t}$ is updated
with all the points $\mathbf{x}$ in $\mathbf{S}_{t-1}$ and $\mathbf{P}_{t-1}$
that have been evaluated during the previous iteration. The process is
summarized in Algorithm 7.

$\begin{array}[]{l}\text{{[1] Initialization}}\\\
\;\;\;\;\begin{array}[]{l}\text{$t\leftarrow 0$}\\\ \text{Set initial poll and
mesh sizes $\Delta^{\text{P}}_{0}\geq\Delta^{\mathcal{M}}_{0}>0$}\\\
\text{Initialize $\mathbf{X}_{0}$ with starting points}\\\ \text{Evaluate
$\\{f(\mathbf{x}),c_{1}(\mathbf{x}),c_{2}(\mathbf{x}),\ldots,c_{m}(\mathbf{x})\\}\;\forall\mathbf{x}\in\mathbf{X}_{0}$}\\\
\end{array}\\\ \text{{[2] Model search}}\\\
\;\;\;\;\begin{array}[]{l}\text{Use $\mathbf{X}_{t}$ to build $\hat{f}$ and
$\\{\hat{c}_{j}\\}_{j\in J}$}\\\ \text{Solve surrogate
problem~{}\eqref{eq:SurrogateProblem} using the inner MADS instance}\\\
\text{$\mathbf{\hat{X}}_{t}\leftarrow$ Set of points evaluated with surrogate
model while solving~{}\eqref{eq:SurrogateProblem}}\\\
\text{$\mathbf{S}_{t}\leftarrow$ Cycle through selection steps to select $q$
points of $\mathbf{\hat{X}}_{t}$}\\\ \text{$\mathbf{S}_{t}\leftarrow$
Projection of the points of $\mathbf{S}_{t}$ onto mesh $\mathcal{M}_{t}$}\\\
\text{Parallel evaluation of
$\\{f(\mathbf{x}),c_{1}(\mathbf{x}),c_{2}(\mathbf{x}),\ldots,c_{m}(\mathbf{x})\\}\;\forall\mathbf{x}\in\mathbf{S}_{t}$}\\\
\text{If success, {\tt goto} {[4]}}\\\ \end{array}\\\ \text{{[3] Poll}}\\\
\;\;\;\;\begin{array}[]{l}\text{Build poll set $\mathbf{P}_{t}$}\\\ \text{Sort
$\mathbf{P}_{t}$ according to $\hat{f}$ and $\\{\hat{c}_{j}\\}_{j\in J}$}\\\
\text{Parallel evaluation of
$\\{f(\mathbf{x}),c_{1}(\mathbf{x}),c_{2}(\mathbf{x}),\ldots,c_{m}(\mathbf{x})\\}\;\forall\mathbf{x}\in\mathbf{P}_{t}$}\\\
\end{array}\\\ \text{{[4] Updates}}\\\
\;\;\;\;\begin{array}[]{l}\text{$t\leftarrow t+1$}\\\ \text{Update
$\Delta^{\mathcal{M}}_{t}$, $\Delta^{\text{P}}_{t}$, $\mathbf{x}^{*}_{t}$ and
$\mathbf{X}_{t}$}\\\ \text{If no stopping condition is met, {\tt goto}
{[2]}}\\\ \end{array}\\\ \end{array}$

Algorithm 7 The proposed MADS optimization algorithm

##  5 Numerical investigation

The proposed SEARCH step technique is tested using five optimization problems.
We first describe the algorithms considered for benchmarking. Next, numerical
results are presented and discussed for the five engineering design problems.

###  5.1 Compared algorithms

Five solvers are compared in our numerical experiments, all based on the MADS
algorithm and implemented using NOMAD 3.8 [3]. This ensures avoiding coding
biases since features are identical among solvers.

  * •

MADS. Refers to the POLL step of MADS, without any SEARCH step, where $2n$
directions are generated and evaluated in parallel. If needed, $k$ additional
directions are generated such that $2n+k$ is a multiple of $q$.

  * •

Multi-Start. Consists of $q$ parallel runs of MADS. They are totally
independent and each instance runs on its own CPU. Each instance proceeds to
its evaluations sequentially, one after the other. Only the POLL step is
executed and no cache is shared between running instances.

  * •

LH Search. The MADS solver mentioned above using a Latin hypercube search
(LHS) at the SEARCH step, where $q$ candidates are generated and evaluated in
parallel.

  * •

Lowess-A. The MADS solver mentioned above with the described surrogate
optimization conducted at the SEARCH step. The $q$ candidates are selected by
cycling through Methods 1 and 2, and then evaluated in parallel.

  * •

Lowess-B. The MADS solver mentioned above with the proposed surrogate
optimization conducted at the SEARCH step. The $q$ search candidates are
selected by cycling through Methods 3, 4, 5, and 6, and then evaluated in
parallel.

Both LOWESS solvers are exactly like Algorithm 7, excepted for the used
selection methods. The only difference between them and the LHS solver is that
the surrogate optimization approach is replaced by a LHS at the SEARCH step.
This should allow us to determine whether surrogate optimization has any
advantage over a random search. The MADS solver is used as the baseline.
Finally, the Multi-Start solver is considered to ensure that one should not
proceed with $q$ independent narrow trajectories instead of one single
trajectory having $q$ wide evaluations.

###  5.2 Engineering design optimization problems

The above solvers are compared on five engineering design application
problems. A short description follows below for each problem. More details are
provided in Appendix B.

  * •

TCSD. The Tension/Compression Spring Design problem consists of minimizing the
weight of a spring under mechanical constraints [24, 25, 26]. This problem has
three variables and four constraints. The design variables define the geometry
of the spring. The constraints concern shear stress, surge frequency, and
minimum deflection.

  * •

Vessel. This problem considers the design of a compressed air storage tank and
has four design variables and four constraints [24, 27]. The variables define
the geometry of the tank and the constraints are related to the volume,
pressure, and solidity of the tank. The objective is to minimize the total
cost of the tank, including material and labour.

  * •

Welded. The welded beam design problem (Version I) has four variables and six
constraints [24, 28]. It aims at minimizing the construction cost of a beam,
under shear stress, bending stress, and deflection constraints. The design
variables define the geometry and the characteristics of the welded joint.

  * •

Solar 1. This optimization problem aims at maximizing the energy received over
a period of 24 hours under five constraints related to budget and heliostat
field area [29]. It has nine variables, including an integer one without an
upper bound.

  * •

Solar 7. This problem aims at maximizing the efficiency of the receiver over a
period of 24 hours for a given heliostats field under six binary constraints
[29]. It has seven design variables, including an integer one without an upper
bound.

A progressive barrier is used to deal with the aggregated constraints [13].
The three first problems are easier relative to the last two ones. However, it
is difficult to find a feasible solution for the TCSD problem. Among all the
considered problems, Solar 1 is certainly the most difficult one.

###  5.3 Numerical experiments

We compare the efficiency of each solver for different values of block size
$q\in\\{1,2,4,8,16,32,64\\}$. As an example, we will use “Lowess-A 16” to
refer to the solver that relies on LOWESS models cycling over Methods 1 and 2
considering a block size $q=16$. For each problem, we generated 50 sets of 64
starting points with Latin hypercube sampling [20]. For all solvers other than
“Multi-Start”, only the first point of each set is used to perform
optimizations. Doing so, we get 50 runs from the same starting points for each
solver, each problem, and each value of $q$. For “Multi-Start”, since $q$
independent and parallel sequential runs of MADS must be performed, we use the
$q$ first points of each set. Doing so, we still get 50 runs for each problem
and each $q$, while ensuring that all starting points are the same for all
solvers.

To avoid that all “LH Search” runs end up with nearly identical solutions for
a given $q$, we use a random seed for initializing each LHS.

For the relatively three simpler problems (TCSD, Vessel, and Welded), a budget
of 100 block evaluations is allocated. For the two relatively difficult
problems (Solar 1 and Solar 7), the budget is increased to 200 block
evaluations. This means that, for a given problem, all solvers will have the
same “wall-clock” time, but not necessarily the same resources (number of CPUs
available for block evaluations) nor the same total number of blackbox
evaluations.

#### Solution quality

Figures 1 and 2 represent the distribution of the final objective function
over the 50 runs for each problem, each solver, and each block size $q$. The
minimum and maximum objective values that we obtained from the runs are
indicated by circles in the figures. Lower and upper quantiles are delimited
by boxes. Median values are represented by a bar into the boxes. The more a
distribution is on the left side, the better the combination of solver and $q$
is. Since we are mostly interested in the best combinations, the figures only
focus on the smallest values. Otherwise, it would be difficult to discern the
difference among the best combinations. All combinations for which the
distribution is cut on the right side are performing poorly.

Figure 1: Performance summary for the TCSD, Vessel, and Welded problems over
50 runs Figure 2: Performance summary for the Solar 1 and Solar 7 problems
over 50 runs

For the three simpler problems (TCSD, Vessel, and Welded, Figure 1), the
LOWESS solvers (and in particular “Lowess-B”) are by far superior to the
solvers that do not rely on surrogate optimization.

For the TCSD problem, the “MADS” solver often failed to find a feasible design
(thus leading to infinite objective values), even with a large number of
evaluations per block. The four other solvers always managed to find a
feasible point for at least 75% of the runs. The “Lowess-B” solver performs
better than any of the other ones. We see that “Lowess-A 64” is outperformed
by “Lowess-B 8”. As the TCSD problem is very constrained, the final objective
function value depends on the initial guess. This is why the “Multi-Start”
solver performs quite well on this problem.

The same trend is observed for the Vessel, and Welded problems. “Lowess-B”
performs better than “Lowess-A”, which outperforms “LHS” or “MADS”. In
particular, “Lowess-B 8” outperforms the solvers “Lowess-A 8/16/32”. As
expected, increasing the block size improves performance. However, for the
“Lowess-B” solver these three problems are easy to solve, so it is difficult
to see an advantage of using parallel computing because the global optimum is
found most of the time within 100 block evaluations for a block size of 16 or
more. The “Multi-Start” solver performs rather poorly on these two problems.

The numerical results generally follow the same trend for the two Solar
problems (Figure 2). For a block of equal size, the LOWESS solvers outperform
the other solvers while “Lowess-B” outperforms “Lowess-A”.

#### Convergence rate

We now examine the convergence rate of the solvers for the case where $q=64$.
Figures 3 and 4 depict the evolution of the median objective function value of
50 runs as a function of the number of block evaluations. For each problem,
the plots on the left compare the convergence of the five solvers with blocks
of size $q=64$ while the plots on the right compare the convergence of the
best-performing solver i.e., “Lowess-B”, for block sizes ranging from $q=1$ to
64.

TCSD problem

Vessel problem

Welded problem

Figure 3: Results for the TCSD, Vessel, and Welded problems; median objective
value of 50 runs

Solar 1

Solar 7

Figure 4: Results for the Solar 1 and Solar 7 problems; median objective value
of 50 runs

We can conclude that “Lowess-B” yields the best solutions faster than any
other solver (for $q=64$). The worst-performing solvers are “Multi-Start” for
problems TCSD and Solar 1 and “LH Search” for problems Vessel, Welded, and
Solar 7. It is also notable that although “MADS” does not use the SEARCH step,
it performs generally well, except for Solar 1. “Lowess-A” performed well but
does not clearly outperform other solvers.

Considering the performance of “Lowess-B” as a function of $q$, we observe
that, as expected, convergence improves for larger values of $q$. Depending on
the problem, there may be a saturation point beyond which an increase of $q$
does not effect an improvement. E.g., a saturation point arises around $q=8$
or 16 for TCSD, Vessel, and Welded. On the contrary, $q$ could be even larger
than 64 for Solar 1 as more CPUs can be utilized.

#### Performance profiles

We now consider performance profiles, which indicate the percentage of runs
where the problem is solved within a deviation from the best known solution
$\tau$ under a budget of function evaluations [30]. Specifically, for each
solver $s$, each instance $r$ and problem $p$, we compute the number of block
evaluations $b_{s,p,r}(\tau)$ such that

$\frac{|f_{s,b,p,r}-f_{p}^{*}|}{|f_{p}^{*}|}\leq\tau,$ (13)

where $f_{p}^{*}$ is the best known objective value for problem $p$ and
$f_{s,b,p,r}$ is the value obtained with the solver $s$ after $b$ block
evaluations. Let $b^{\min}_{p,r}(\tau)$ be the smallest budget for solving the
instance $r$ of problem $p$ with deviation $\tau$, i.e.,

$b^{\min}_{p,r}(\tau)=\underset{s}{\min}\;b_{s,p,r}(\tau).$

Then, we can plot the proportion of runs of solver $s$ that satisfy Eq. (13)
at a multiple $\alpha$ of the smallest budget, i.e.,
$\alpha\,b^{\min}_{p,r}(\tau)$ block evaluations. Figure 5 depicts the
performance profiles of the five considered problems for $q=64$ over the 50
runs (instances) for $\tau$ values that range from 10-1 to 10-4.

Performance profiles for $\tau=10^{-1}$ | Performance profiles for $\tau=10^{-2}$
---|---
|
Performance profiles for $\tau=10^{-3}$ | Performance profiles for $\tau=10^{-4}$
|
Figure 5: Performance profiles for $q=64$ over the 50 runs of all five
problems

Note that higher curves imply better solver performance. Moreover, a solver
that performs well for small values of $\alpha$ is a solver that can solve,
for the considered $\tau$, a large number of problems with a small evaluation
budget. Figure 5 confirms our previous observations: “Lowess-B” outperforms
all solvers, followed by “Lowess-A” and “LH Search” in second and third
position. “MADS” and “Multi-Start” are in the last position. For large
tolerances, “MADS” does better than “Multi-Start”, and better than “LH Search”
for small values of $\alpha$ ($\leq 4$). For small precisions, “MADS” is
outperformed by the other solvers.

The most interesting observation is the significant gap between the
performance curves of LOWESS solvers and the ones from the three other
solvers. The gap increases as $\tau$ decreases. From moderate to lower $\tau$
values ($\leq 10^{-2}$), “Lowess-B” systematically solves at least twice more
problems than “LH Search”, “Multi-Start” and “MADS”. It is unusual to observe
such clear differences on performance profiles. “Lowess-A” performs almost as
well as “Lowess-B”, particularly when $\tau$ is small (around $10^{-4}$), but
needs at least four times more block evaluations to achieve this.

#### Scalability analysis

We wish to establish the reduction of wall-clock time when using additional
resources for each solver. To that end, we follow the methodology proposed in
[31]. We define $f_{s,b,p,r,q}$ as the value of the objective function
obtained by solver $s$ after $b$ block evaluations on instance $r$ of problem
$p$ when using $q$ CPUs. We also define the _reference_ objective value as the
best value achieved with only one CPU ($q=1$), i.e.,

$f^{\text{ref}}_{s,p,r}=\underset{b}{\min}f_{s,b,p,r,1},$ (14)

and $b^{\text{ref}}_{s,p,r,q}$ the number of block evaluations necessary to
reach $f^{\text{ref}}_{s,p,r}$ when $q$ CPUs are used, i.e.,

$b^{\text{ref}}_{s,p,r,q}=\min\\{b:f^{\text{ref}}_{s,p,r}\leq
f_{s,b,p,r,q}\\}.$ (15)

The _speed-up_ of solver $s$ when solving with $q$ CPUs is defined as

$\text{speed-
up}(s,q)=\underset{p,r}{\text{geomean}}\left(\;\frac{b^{\text{ref}}_{s,p,r,q}}{b^{\text{ref}}_{s,p,r,1}}\;\right)$
(16)

and its _efficiency_ as

$\text{efficiency}(s,q)=\frac{\text{speed-up}(s,q)}{q}.$ (17)

Figure 6 depicts the speed-up and efficiency values obtained by our numerical
experiments.

Figure 6: Speed-up and efficiency

Perfect scalability is obtained when the speed-up is equal to $q$ and the
efficiency is equal to 1. The speed-up curves show that the power introduced
by new CPUs decreases as their number increases. This was observed in Figures
3 and 4 where problems exhibited saturation around $q=8$ and 16. “Lowess-B”
achieves the best speed-up, followed by “MADS”. “Lowess-A” and “LH Search”
come next, followed by the worst-performing solver, namely “Multi-Start”. We
conclude that it is better and more productive to proceed with one search
performing $q$ parallel evalutions instead of conducting $q$ independent
searches consisting of a single evaluation.

The efficiency curves demonstrate rapid decrease except for “Lowess-B”; its
rate exhibits a bump on its efficiency curve at $q=4$. For $q=2$, only methods
3 and 4 are used to generate candidates. For $q\geq 4$, methods 5 and 6 are
also used. The aforementioned bump highlights the important contributions of
these methods to the efficiency of “Lowess-B”.

##  6 Conclusion

Linear LOWESS models with optimized kernel shapes and coefficients seem to
provide high-performing surrogates of the blackboxes. The use of diverse
selection methods (3 to 6) enables an efficient exploration of the design
space, accelerates local convergence, and makes optimal use of additional CPU
resources. Methods 5 and 6 are particularly efficient, outperforming the other
selection methods. This means that the way surrogates are used by method 1 is
not effective. Similarly, the diversification strategy of method 2 is not
adequate to select points that lie far enough from the ones already evaluated.

We cannot draw a definite conclusion about which of the methods 5 or 6 is
better than the other. We believe that the good performance of “Lowess-B” is
due to using method 5.

The proposed selection methods are not specific to the LOWESS model considered
here; they are applicable to any surrogates. We believe that they will work
well with reduced-fidelity (or variable-fidelity) physical-based models since
high- and low-fidelity models typically have similarly structured solution
domains. The selection methods are also applicable to other algorithms using
surrogates to identify promising points to evaluate.

## References

  * [1] O. Kramer, D.E. Ciaurri, and S. Koziel.  Derivative-free optimization.  In S. Koziel and XS. Yang, editors, Computational Optimization, Methods and Algorithms, volume 356 of Studies in Computational Intelligence, pages 61–83. Springer, 2011.
  * [2] C. Audet, S. Le Digabel, C. Tribes, and V. Rochon Montplaisir.  The NOMAD project.  Software available at https://www.gerad.ca/nomad.
  * [3] S. Le Digabel.  Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm.  ACM Transactions on Mathematical Software, 37(4):44:1–44:15, 2011\.
  * [4] C. Audet and J.E. Dennis, Jr.  Mesh adaptive direct search algorithms for constrained optimization.  SIAM Journal on Optimization, 17(1):188–217, 2006.
  * [5] E. Alba, G. Luque, and S. Nesmachnow.  Parallel metaheuristics: recent advances and new trends.  International Transactions in Operational Research, 20(1):1–48, 2013.
  * [6] S. Le Digabel, M.A. Abramson, C. Audet, and J.E. Dennis, Jr.  Parallel versions of the MADS algorithm for black-box optimization.  In Optimization days, Montreal, May 2010. GERAD.  Slides available at http://www.gerad.ca/Sebastien.Le.Digabel/talks/2010_JOPT_25mins.pdf.
  * [7] B. Talgorn, S. Le Digabel, and M. Kokkolaras.  Statistical Surrogate Formulations for Simulation-Based Design Optimization.  Journal of Mechanical Design, 137(2):021405–1–021405–18, 2015\.
  * [8] Mahdi Pourbagian, Bastien Talgorn, WagdiG. Habashi, Michael Kokkolaras, and Sébastien Le Digabel.  Constrained problem formulations for power optimization of aircraft electro-thermal anti-icing systems.  Optimization and Engineering, pages 1–31, 2015.
  * [9] C. Audet, M. Kokkolaras, S. Le Digabel, and B. Talgorn.  Order-based error for managing ensembles of surrogates in derivative-free optimization.  Journal of Global Optimization, 70(3):645–675, 2018.
  * [10] B. Talgorn, C. Audet, M. Kokkolaras, and S. Le Digabel.  Locally weighted regression models for surrogate-assisted design optimization.  Optimization and Engineering, 19(1):213–238, 2018.
  * [11] Raphael T. Haftka, Diane Villanueva, and Anirban Chaudhuri.  Parallel surrogate-assisted global optimization with expensive functions – a survey.  Structural and Multidisciplinary Optimization, 54(1):3–13, Jul 2016\.
  * [12] R. Fletcher and S. Leyffer.  Nonlinear programming without a penalty function.  Mathematical Programming, Series A, 91:239–269, 2002.
  * [13] C. Audet and J.E. Dennis, Jr.  A progressive barrier for derivative-free nonlinear programming.  SIAM Journal on Optimization, 20(1):445–472, 2009.
  * [14] R. Alizadeh, J.K. Allen, and F. Mistree.  Managing computational complexity using surrogate models: a critical review.  Research in Engineering Design, 2020.
  * [15] W.S. Cleveland.  Robust locally weighted regression and smoothing scatterplots.  Journal of the American Statistical Association, 74:829–836, 1979\.
  * [16] W.S. Cleveland.  LOWESS: A Program for Smoothing Scatterplots by Robust Locally Weighted Regression.  The American Statistician, 35(1), 1981.
  * [17] W.S. Cleveland and S.J. Devlin.  Locally weighted regression: An approach to regression analysis by local fitting.  Journal of the American Statistical Association, 83:596–610, 1988\.
  * [18] W.S. Cleveland, S.J. Devlin, and E. Grosse.  Regression by local fitting: methods, properties, and computational algorithms.  Journal of Econometrics, 37(1):87 – 114, 1988.
  * [19] B. Talgorn.  SGTELIB: Surrogate model library for derivative-free optimization.  https://github.com/bbopt/sgtelib, 2019.
  * [20] M.D. McKay, R.J. Beckman, and W.J. Conover.  A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.  Technometrics, 21(2):239–245, 1979.
  * [21] T.J. Santner, B.J. Williams, and W.I. Notz.  The Design and Analysis of Computer Experiments, chapter 5.2.2, Designs Generated by Latin Hypercube Sampling, pages 127–132.  Springer, New York, NY, 2003.
  * [22] N. Mladenović and P. Hansen.  Variable neighborhood search.  Computers and Operations Research, 24(11):1097–1100, 1997.
  * [23] P. Hansen and N. Mladenović.  Variable neighborhood search: principles and applications.  European Journal of Operational Research, 130(3):449–467, 2001\.
  * [24] H. Garg.  Solving structural engineering design optimization problems using an artificial bee colony algorithm.  Journal of Industrial and Management Optimization, 10(3):777–794, 2014.
  * [25] J. Arora.  Introduction to Optimum Design.  Elsevier Science, 2004.
  * [26] A.D. Belegundu.  A Study of Mathematical Programming Methods for Structural Optimization.  University of Iowa, 1982.
  * [27] B. K. Kannan and S. N. Kramer.  Augmented Lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design.  Journal of Mechanical Design, 65:103–112+, 1993.
  * [28] Singiresu S. Rao.  Engineering Optimization: Theory and Practice, 3rd Edition.  Wiley-Interscience, 1996.
  * [29] Mathieu Lemyre Garneau.  Modelling of a solar thermal power plant for benchmarking blackbox optimization solvers.  Master’s thesis, École Polytechnique de Montréal, 2015.
  * [30] E.D. Dolan and J.J. Moré.  Benchmarking optimization software with performance profiles.  Mathematical Programming, 91(2):201–213, 2002.
  * [31] Prasad Jogalekar and Murray Woodside.  Evaluating the scalability of distributed systems.  IEEE Trans. Parallel Distrib. Syst., 11(6):589–603, June 2000\.
  * [32] C.G. Atkeson, A.W. Moore, and S. Schaal.  Locally weighted learning.  Artificial Intelligence Review, pages 11–73, 1997.

##  Appendix A LOWESS predictions

As a convention, we denote with
${\boldsymbol{\xi}}\in\mathcal{X}\subseteq{\mathbb{R}}^{n}$ the point of the
design space where we want to predict the value of the blackbox output.
Locally weighted scatterplot smoothing (LOWESS) models build a local linear
regression at the point ${\boldsymbol{\xi}}$ where the blackbox output
$[f\;c_{1}\ldots c_{m}]$ are to be estimated [32, 15, 16, 17, 18, 10]. This
local regression emphasizes data points that are close to
${\boldsymbol{\xi}}$. The interested reader can refer to [10] for details
about the method described below. We consider here only local linear
regressions; local quadratic regressions and Tikhonov regularization are
considered in [10]. On the contrary, while only a Gaussian kernel was
considered in [10], six others are added here as kernel functions.

We define the output matrix $\mathbf{Y}\in{\mathbb{R}}^{p\times(m+1)}$, the
design matrix $\mathbf{Z}_{\boldsymbol{\xi}}\in{\mathbb{R}}^{p\times(n+1)}$,
and the weight matrix $\mathbf{W}_{\boldsymbol{\xi}}\in{\mathbb{R}}^{p\times
p}$:

$\mathbf{Y}=\left[\begin{array}[]{c c @{\,} c @{\,}
c}f(\mathbf{x}_{1})&c_{1}(\mathbf{x}_{1})&\ldots&c_{m}(\mathbf{x}_{1})\\\
\vdots&\vdots&&\vdots\\\
f(\mathbf{x}_{p})&c_{1}(\mathbf{x}_{p})&\ldots&c_{m}(\mathbf{x}_{p})\end{array}\right],\;\;\mathbf{Z}_{\boldsymbol{\xi}}=\left[\begin{array}[]{c
c}1&(\mathbf{x}_{1}-{\boldsymbol{\xi}})^{\top}\\\ \vdots&\vdots\\\
1&(\mathbf{x}_{p}-{\boldsymbol{\xi}})^{\top}\end{array}\right],\;\;\mathbf{W}_{\boldsymbol{\xi}}=\left[\begin{array}[]{c
@{} c @{} c}w_{1}({\boldsymbol{\xi}})&&\\\ &\ddots&\\\
&&w_{p}({\boldsymbol{\xi}})^{\top}\end{array}\right].$ (18)

The details of the computation of $w_{i}({\boldsymbol{\xi}})$ are described in
Section A.1. Then, we define
$\mathbf{u}_{\boldsymbol{\xi}}\in{\mathbb{R}}^{n+1}$ as the first column of
$(\mathbf{Z}_{\boldsymbol{\xi}}^{\top}\mathbf{W}_{\boldsymbol{\xi}}\mathbf{Z}_{\boldsymbol{\xi}})^{-1}$,
which means that $\mathbf{u}_{\boldsymbol{\xi}}$ is the solution of the linear
system
$\mathbf{Z}_{\boldsymbol{\xi}}^{\top}\mathbf{W}_{\boldsymbol{\xi}}\mathbf{Z}_{\boldsymbol{\xi}}\mathbf{u}_{\boldsymbol{\xi}}=\mathbf{e}_{1}$.
The prediction of the blackbox outputs at ${\boldsymbol{\xi}}$ is then

$\hat{\mathbf{y}}({\boldsymbol{\xi}})=\left[\begin{array}[]{c c c
c}\hat{f}({\boldsymbol{\xi}})&\hat{c}_{1}({\boldsymbol{\xi}})&\ldots&\hat{c}_{m}({\boldsymbol{\xi}})\end{array}\right]=\mathbf{u}_{\boldsymbol{\xi}}^{\top}\mathbf{Z}_{\boldsymbol{\xi}}^{\top}\mathbf{W}_{\boldsymbol{\xi}}\mathbf{Y}.$
(19)

The cross-validation value
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{y}}}(\mathbf{x}_{i})$
(i.e., the value of the LOWESS model at $\mathbf{x}_{i}$ when the data point
$\mathbf{x}_{i}$ is not used to build the model) are computed by setting
$w_{i}$ to 0. Unfortunately, unlike for radial basis function (RBF) models or
polynomial response surfaces (PRSs), we do not know any computational shortcut
allowing a more efficient computation of the values of
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{y}}}$. However, each value
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{y}}}(\mathbf{x}_{i})$ is
computed at the same computational cost as a prediction
$\hat{y}(\mathbf{x}_{i})$.

###  A.1 Weights computation in LOWESS models

The weight $w_{i}({\boldsymbol{\xi}})$ quantifies the relative importance of
the data point $\mathbf{x}_{i}$ in the construction of the local regression at
${\boldsymbol{\xi}}$. Like for kernel smoothing, it relies on a kernel
function $\phi$ and depends on the distance between ${\boldsymbol{\xi}}$ and
$\mathbf{x}_{i}$. In our method, we use

$w_{i}({\boldsymbol{\xi}})=\phi\left(\lambda\frac{\|{\boldsymbol{\xi}}-\mathbf{x}_{i}\|_{2}}{d_{n+1}({\boldsymbol{\xi}})}\right),$
(20)

where $\phi(d)$ is one of the kernel functions described in Table 1 and Figure
7. All kernel functions are normalized so that $\phi(0)=1$ and, if applicable,
$\int_{\mathbb{R}}\phi=1$. As the integral of the inverse multi-quadratic
kernel does not converge, the normalization constant 52.015 is introduced to
minimize the $\mathcal{L}^{2}$ distance between the inverse multi-quadratic
and inverse quadratic kernel. The parameter $\lambda>0$ controls the general
shape of the model, and $d_{n+1}({\boldsymbol{\xi}})$ is a local scaling
coefficient that estimates the distance of the $n+1^{th}$ closest data point
to ${\boldsymbol{\xi}}$. The kernel function $\phi$ and the shape parameter
$\lambda$ are chosen to minimize the aggregate order error with cross-
validation (AOECV) described in Section A.2. The fact that some of the
available kernel function have a compact domain gives to LOWESS models the
ability to ignore outliers or aberrant data points. As an example, if the
blackbox fails to compute correctly the objective function for a given data
point, the value returned by the blackbox might be an arbitrarily high value
(e.g., $1.8~{}10^{308}$ for a C++ code returning the standard max double).
With non-compact kernel function, this would perturb the LOWESS model on the
entire design space. However, if there is no such aberrant data points, non-
compact kernel functions tend to yield better results.

Table 1: Possible values for the kernel function $\phi$ # | Kernel name | $\phi:{\mathbb{R}}\rightarrow{\mathbb{R}}^{+}$ | Compact domain
---|---|---|---
1 | Tri-cubic | $\phi(d)=(1-|\frac{162}{140}d|^{3})^{3}\mathbb{1}_{|d|\leq\frac{140}{162}}$ | Yes
2 | Epanechnikov | $\phi(d)=(1-\frac{16}{9}d^{2})\mathbb{1}_{|d|\leq\frac{3}{4}}$ | Yes
3 | Bi-quadratic | $\phi(d)=(1-|\frac{16}{15}d|^{2})^{2}\mathbb{1}_{|d|\leq\frac{15}{16}}$ | Yes
4 | Gaussian | $\phi(d)=\exp(-\pi d^{2})$ | No
5 | Inverse quadratic | $\phi(d)=\frac{1}{1+\pi^{2}d^{2}}$ | No
6 | Inverse multi-quadratic | $\phi(d)=\frac{1}{\sqrt{1+52.015d^{2}}}$ | No
7 | Exp-root | $\phi(d)=\exp(-2\sqrt{|d|})$ | No
Figure 7: Representation of the 7 kernels listed in Table 1

To obtain a model $\hat{y}$ that is differentiable everywhere, [10] defines
$d_{n+1}({\boldsymbol{\xi}})$ such that the expected number of training points
in a ball of center ${\boldsymbol{\xi}}$ and radius
$d_{n+1}({\boldsymbol{\xi}})$ is $n+1$:

${\mathbb{E}}\left[\text{card}\Big{\\{}\mathbf{x}_{i}:\mathbf{x}_{i}\in\mathbf{X},\|{\boldsymbol{\xi}}-\mathbf{x}_{i}\|_{2}\leq
d_{n+1}({\boldsymbol{\xi}})\Big{\\}}\right]=n+1.$ (21)

Moreover, [10] observes that the values
$\big{\\{}\|{\boldsymbol{\xi}}-\mathbf{x}_{i}\|_{2}^{2}\big{\\}}_{i=1,\ldots,p}$
can be fitted well by a Gamma distribution and therefore defines the local
scaling parameter as

$d_{n+1}({\boldsymbol{\xi}})=\sqrt{g^{(-1)}\left(\frac{\mu_{\boldsymbol{\xi}}^{2}}{\sigma_{\boldsymbol{\xi}}^{2}},\frac{\sigma_{\boldsymbol{\xi}}^{2}}{\mu_{\boldsymbol{\xi}}};\frac{n+1}{p}\right)},$
(22)

where $\mu_{\boldsymbol{\xi}}$ (resp. $\sigma^{2}_{\boldsymbol{\xi}}$) denotes
the mean (resp. variance) of $\|{\boldsymbol{\xi}}-\mathbf{x}_{i}\|_{2}^{2}$
over $\mathbf{X}$ and $g^{(-1)}(k,\theta;.)$ is the inverse function of the
cumulative density function of a Gamma distribution with shape parameter $k$
and scale parameter $\theta$.

###  A.2 Aggregate Order Error with Cross-Validation

The AOECV is an error metric that aims at quantifying the quality of a _multi-
output_ surrogate model. Specifically, it aims at quantifying the discrepancy
between problems ($P$) and ($\hat{P}$) for a given surrogate model. We first
define the aggregate constraint violation function [12]
$h(\mathbf{x})=\sum_{j=1}^{m}\max\\{0,c_{j}(\mathbf{x})\\}^{2}$. Note that
other definitions of $h$ are possible (notably: number of violated
constraints, most violated constraint, etc.) but as the previous definition of
$h$ is used in the main MADS instance to solve ($P$), we need to use the same
aggregate constraint in our definition of the AOECV.

We then define the order operators

$\displaystyle\mathbf{x}\prec\mathbf{x}^{\prime}\Leftrightarrow$
$\displaystyle\left\\{\begin{array}[]{l}h(\mathbf{x})<h(\mathbf{x}^{\prime})\\\
\text{or}\\\ h(\mathbf{x})=h(\mathbf{x}^{\prime})\text{ and
}f(\mathbf{x})<f(\mathbf{x}^{\prime}),\end{array}\right.$ (26)
$\displaystyle\mathbf{x}\preceq\mathbf{x}^{\prime}\Leftrightarrow$
$\displaystyle\;\;{\tt not}(\mathbf{x}^{\prime}\prec\mathbf{x})$ (27)

which are transitive. In particular, the incumbent solution $\mathbf{x}^{t}$
of the original problem ($P$) is such that
$\mathbf{x}^{t}\preceq\mathbf{x},\;\forall\mathbf{x}\in\mathbf{X}$. Similarly,
a global minimizer $\mathbf{x}^{*}$ is such that
$\mathbf{x}^{*}\preceq\mathbf{x},\;\forall\mathbf{x}\in\mathcal{X}$. By the
same principle, we define the operator
$~{}\widehat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\widehat{\prec}}}~{}$ by
using the cross-validation values
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{f}}}$ and
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{h}}}=\sum_{j=1}^{m}\max\\{0,\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{c}}}_{j}(\mathbf{x})\\}^{2}$
instead of $f$ and $h$. We then define the aggregated order error with cross-
validation (AOECV) metric:

$\mathcal{E}_{AOECV}=\frac{1}{p^{2}}\displaystyle\sum_{i=1}^{p}\sum_{j=1}^{p}~{}{\tt
xor}~{}\Big{(}\mathbf{x}_{i}\prec\mathbf{x}_{j},\mathbf{x}_{i}~{}\widehat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\widehat{\prec}}}~{}\mathbf{x}_{j}\Big{)}.$
(28)

where xor is the exclusive or operator (i.e., $~{}{\tt xor}~{}(A,B)=1$ if the
booleans $A$ and $B$ differ and 0 otherwise). The metric allows to quantify
how often the model is able to correctly decide which of two points is better.

The shape parameter $\lambda$ and the kernel function $\phi$ are then chosen
to minimize $\mathcal{E}_{AOECV}(\lambda,\phi).$ If two couples
$(\lambda,\phi)$ lead to the same metric value (because of the piecewise-
constant nature of the metric), the couple with the smallest value of
$\lambda$ (i.e., the smoother model) is preferred.

##  Appendix B Detailed description of the test problems

The five engineering design application problems considered are listed in
Table 2. Problem size is reflected by $n$ and $m$, where $n$ denotes the
number of design variables and $m$ the number of general nonlinear inequality
constraints. Table 2 also indicates whether any variables are integer or
unbounded and reports the best known value of the objective function.

Table 2: Summary of the five engineering design optimization problems Problem | $n$ | $m$ | Integer | Infinite | Best objective
---|---|---|---|---|---
name |  |  | variables | bounds | function value
TCSD | 3 | 4 | No | No | 0.0126652
Vessel | 4 | 4 | No | No | 5,885.332
Welded | 4 | 6 | No | No | 2.38096
Solar 1 | 9 | 5 | Yes | Yes | –900,417
Solar 7 | 7 | 6 | Yes | Yes | –4,976.17

The Tension/Compression Spring Design (TCSD) problem consists of minimizing
the weight of a spring under mechanical constraints [24, 25, 26]. The design
variables define the geometry of the spring. The constraints concern shear
stress, surge frequency and minimum deflection. The best known solution,
denoted $\mathbf{x}^{*}$, and the bounds on the variables, denoted by
$\underline{\mathbf{x}}$ and $\bar{\mathbf{x}}$, are given in Table 3.

Table 3: Variables of the TCSD problem Variable description | $\underline{\mathbf{x}}$ | $\bar{\mathbf{x}}$ | $\mathbf{x}^{*}$
---|---|---|---
Mean coil diameter | 0.05 | 2 | 0.051686696913218
Wire diameter | 0.25 | 1.3 | 0.356660815351066
Number of active coil | 2 | 15 | 11.292312882259289

The Vessel problem considers the optimal design of a compressed air storage
tank [24, 27]. The design variables define the geometry of the tank The
constraints are related to the volume, pressure, and solidity of the tank. The
objective is to minimize the total cost of the tank, including material and
labour. Table 4 lists the variable bounds and the best known solution.

Table 4: Variables of the Vessel problem Variable description | $\underline{\mathbf{x}}$ | $\bar{\mathbf{x}}$ | $\mathbf{x}^{*}$
---|---|---|---
Thickness of the vessel | 0.0625 | 6.1875 | 0.778168641330718
Thickness of the head | 0.0625 | 6.1875 | 0.384649162605973
Inner radius | 10 | 200 | 40.319618721803231
Length of the vessel without heads | 10 | 200 | 199.999999998822659

The Welded (or welded beam design) problem (Version I) consists of minimizing
the construction cost of a beam under shear stress, bending stress, load and
deflection constraints [24, 28]. The design variables define the geometry of
the beam and the characteristics of the welded joint. Table 5 lists the
variable bounds and the best known solution.

Table 5: Variables of the Welded problem Variable description | $\underline{\mathbf{x}}$ | $\bar{\mathbf{x}}$ | $\mathbf{x}^{*}$
---|---|---|---
Thickness of the weld | 0.1 | 2 | 0.244368407428265
Length of the welded joint | 0.1 | 10 | 6.217496713101864
Width of the beam | 0.1 | 10 | 8.291517255567012
Thickness of the beam | 0.1 | 2 | 0.244368666449562

The Solar1 and Solar7 problems consider the optimization of a solar farm,
including the heliostat field and/or the receiver [29]. The Solar1
optimization problem aims at maximizing the energy received over a period of
24 hours under several constraints of budget and heliostat field area. This
problem has one integer variable that has no upper bound. Table 6 lists the
variable bounds and the best known solution.

Table 6: Variables of the Solar1 problem Variable description | $\underline{\mathbf{x}}$ | $\bar{\mathbf{x}}$ | $\mathbf{x}^{*}$
---|---|---|---
Heliostat height | 1 | 40 | 6.165258994385601
Heliostat width | 1 | 40 | 10.571794049143792
Tower height | 20 | 250 | 91.948461670428486
Receiver aperture height | 1 | 30 | 6.056202026704944
Receiver aperture width | 1 | 30 | 11.674984434929991
Max number of heliostats (Integer) | 1 | $+\infty$ | 1507
Field maximum angular span | 1 | 89 | 51.762281627953051
Minimum distance to tower | 0.5 | 20 | 1.347318830713629
Maximum distance to tower | 1 | 20 | 14.876940809562798

The Solar7 problem aims at maximizing the efficiency of the receiver over a
period of 24 hours, for a given heliostats field, under 6 binary constraints
[29]. This problem has one integer variable that has no upper bound. The
objective function is the energy transferred to the molten salt. Table 7 lists
the variable bounds and the best known solution.

Table 7: Variables of the Solar7 problem Variable description | $\underline{\mathbf{x}}$ | $\bar{\mathbf{x}}$ | $\mathbf{x}^{*}$
---|---|---|---
Aperture height | 1 | 30 | 11.543687848308958
Aperture width | 1 | 30 | 15.244236061098078
Outlet temperature | 793 | 995 | 803.000346734710888
Number of tubes (Integer) | 1 | $+\infty$ | 1292
Insulation thickness | 0.01 | 5 | 3.399190219909724
Tubes inside diameter | 0.005 | 0.1 | 0.010657067457678
Tubes outside diameter | 0.0055 | 0.1 | 0.011167646941518