[1]organization=Lawrence Livermore National Laboratory,country=USA
[2]organization=Department of Mechanical Sciences and Engineering, University
of Illinois at Urbana-Champaign, country=USA

# Design optimization in unilateral contact using pressure constraints and
Bayesian optimization

Jingyi Wang  Jerome Solberg  Mike A. Puso  Eric B. Chin  Cosmin G. Petra

###### Abstract

Design optimization problems, e.g., shape optimization, that involve
deformable bodies in unilateral contact are challenging as they require robust
contact solvers, complex optimization methods that are typically gradient-
based, and sensitivity derivations. Notably, the problems are nonsmooth,
adding significant difficulty to the optimization process. We study design
optimization problems in frictionless unilateral contact subject to pressure
constraints, using both gradient-based and gradient-free optimization methods,
namely Bayesian optimization. The contact simulation problem is solved via the
mortar contact and finite element methods. For the gradient-based method, we
use the direct differentiation method to compute the sensitivities of the cost
and constraint function with respect to the design variables. Then, we use
Ipopt to solve the optimization problems. For the gradient-free approach, we
use a constrained Bayesian optimization algorithm based on the standard
Gaussian Process surrogate model. We present numerical examples that control
the contact pressure, inspired by real-life engineering applications, to
demonstrate the effectiveness, strengths and shortcomings of both methods. Our
results suggest that both optimization methods perform reasonably well for
these nonsmooth problems.

##  1 Introduction

Contacting bodies are ubiquitous in engineering problems, covering, e.g.,
crash worthiness test, part assembly, metal forming [14], etc. Many of these
problems are solved with finite element methods [32]. However, combining
design optimization and finite element methods to _design_ contacting systems
rather than _analyze_ them received little attention. This is partially due to
the complexity of combining state-of-the-art mortar contact solvers with
optimization algorithms and the nonsmoothness that exists in such problems
[11].

Finite element contact solvers have evolved over the years and many techniques
have been developed to solve contact problems in a robust and efficient manner
[10, 22]. Among them, the combined Uzawa-augmented Lagrangian method (ALM) has
proven reliable. Recently, researchers have used the interior-point method, to
solve frictionless contact problems [28]. The constraint function, i.e., gap
function, can be computed using node-on-node, node-on-segment or mortar
method. The mortar segment-on-segment method is preferred because it passes
the patch test and eliminates element locking [22, 23, 35]. In this work, we
use the popular interior-point library Ipopt and the mortar method to solve
the contact problem.

Previous shape and topology optimization work has applied gradient-based
optimization methods to design systems with contact. In [7], the authors
applied topology optimization to design nonlinear elastic bodies in contact
using Uzawa-ALM and mortar methods, adjoint sensitivities and Ipopt. In [6],
the authors used the ALM method to deal with stress-based constraints and
contact boundary conditions. The nonsmoothness of the solution to the contact
problems is highlighted as the reason for the design to be highly dependent on
initial design and optimization algorithms. In [20], the design problems of
sheets in contact are shown to be nonsmooth and solved approximately through a
subgradient method. In [2], a stress-based approach is presented for the
design of truss-like elastic structures through imposing stress constraints.
Sensitivities and gradient-based methods are used in topology optimization
problems in unilateral contact in [26]. In [19], the authors study
crashworthiness optimization of frame structures using a topology optimization
formulation and sensitivities. Contact pressure distribution in plane stress
with friction is optimized in [13] through the smoothing of a pressure
function. While the study is limited to obstacle contact problems, the authors
argued for the necessity of compliance constraint in addition to the objective
function consisting of the mean and variance of pressure. In [17], the authors
solved shape optimization with uniform contact pressure constraints. While the
potential contact surface is not allowed to change, the actual contact region
is affected by the shape design. Though nonsmoothness caused by changing
contact region is not discussed in the paper, the results show promise for
such problems in two dimension with method of moving asymptotes [27]. It is
worth pointing out that most papers on design optimization in contact use
node-on-segment contact or other simplified treatment of contact analysis, and
not mortar contact formulation, rendering the methods less accurate. Moreover,
design objective or constraints that are explicit in pressure are rare.

Despite the success and advancement of these important works, there remain a
few challenges that hinder the further development of a robust and efficient
gradient-based optimizer. The first one is the problem-dependent sensitivity
computation based on complex contact solvers. Not only are derivation and
implementation of sensitivities time-consuming, the necessary information to
compute them is often difficult to obtain. For instance, many design variables
change the shape and geometry of the bodies and thus the discretized mesh. The
sensitivities of the mesh with respect to the design variables might be
unavailable, particularly if third-party mesh tools are employed. Second, the
design optimization problem in contact can be nonsmooth, especially when
considering explicit pressure objective and constraints [11], leaving
researchers with little convergence guarantee even with the most successful
optimizers such as Ipopt or MMA [36]. Here, a smooth function means that it is
continuously differentiable with Lipschitz first-order derivatives [18].
During the optimization, the surface region in contact constantly changes with
updated design variables, which could render the objective and constraint
functions non-differentiable.

To address the first challenge, we propose using the Bayesian optimization
method for design problems in contact. Bayesian optimization is a gradient-
free method that treats the objective and constraints as ‘black-box’ functions
that can be approximated with surrogate models, e.g., Gaussian Process (GP)
[8]. Bayesian optimization has been widely applied to structural design [16],
process optimization of additive manufacturing [30], inertial confinement
fusion design [29], machine learning [1], etc. In [12], a gradient-free
subproblem approximation method is used for the optimization of an insole for
lower plantar fascia stress. To build the surrogate model, at each iteration,
a sample point is collected by running the finite element simulation, where no
sensitivity computation is necessary. Yet Bayesian optimization, particularly
in the constrained case, has its own shortcomings. Without sensitivities, the
algorithm is less likely to produce a local minimum that satisfies first-order
optimality conditions [18]. In addition, for large scale problems, while the
sampling can be easily parallelized, the surrogate model, e.g., Gaussian
Processes (GPs) can become computationally expensive to train with increasing
number of samples.

To study the second challenge, we derive, design and implement sensitivities
for problems with pressure objective and constraints, and solve them using the
well-established interior-point solver Ipopt. The sensitivities are computed
based on the differentiation of the Karush–Kuhn–Tucker (KKT) systems of the
contact solutions. While no new algorithm is proposed to provide convergence
theory, we demonstrate the validity of the gradient-based approach through
numerical examples with changing contact region and pressure constraints that
are inspired by real-life applications. Further, we compare both gradient-
based and gradient-free approaches on our numerical examples. Our main
contributions are

  * 1.

derivation and implementation of a gradient-based method using sensitivities
for design optimization in unilateral contact with pressure constraints

  * 2.

Design and implementation of constrained Bayesian optimization on design
optimization in unilateral contact

  * 3.

Novel design and insights produced by our approach in numerical examples.

The paper is organized as follows. In section 2, we briefly discuss the finite
element contact mechanics models in the optimization form. In section 3, we
formulate the design optimization problem and its sensitivities based on the
finite element model. We propose the constrained Bayesian optimization method
in section 4. The numerical examples are presented in section 5. Finally, the
conclusion is given in section 6.

##  2 Contact problem in optimization formulation

By way of introduction, consider a solid body $B$ in the Euclidean space
$\mathbb{E}^{n}$ and let ${\bm{X}}$ denote the location of its material points
in the undeformed/reference configuration $\Omega$ with boundary
$\partial\Omega$. The body is subjected to an external load field $f_{ext}$
and as such the same material points are displaced to locations in the
deformed configuration $\omega$ such that
${\bm{x}}={\bm{X}}+{\bm{u}}({\bm{X}})$, where
${\bm{u}}:\Omega\to\mathbb{E}^{3}$. The deformation gradient tensor ${\bm{F}}$
is defined as ${\bm{F}}=\mathbf{1}+\frac{\partial{\bm{u}}}{\partial{\bm{X}}}$,
where $\mathbf{1}$ is the mixed identity tensor. We assume $B$ undergoes
infinitesimal deformations and is comprised of an elastic material with
internal energy function $\phi_{e}({\bm{F}},{\bm{x}})$. We also assume the
external force $f_{ext}$ is conservative such that it is derivable from an
external energy function $\phi_{ext}({\bm{u}},{\bm{X}})$.

Now consider two distinct bodies $\Omega^{(1)}$ and $\Omega^{(2)}$ with
material points ${\bm{X}}^{(i)},i=1,2$ that are potentially in contact. The
internal and external energies for $\Omega^{(i)}$ are
$\phi^{(i)}_{e}({\bm{F}}^{(i)},{\bm{X}}^{(i)})$ and
$\phi_{ext}^{(i)}({\bm{u}}^{(i)},{\bm{X}}^{(i)})$.

The boundary $\partial\Omega^{(i)}$ of $\Omega^{(i)}$ is decomposed into the
complementary surfaces $\Gamma_{u}^{i},\Gamma_{t}^{i}$ and $\Gamma_{c}^{i}$ of
Figure 1. Displacement is prescribed over $\Gamma_{u}^{(i)}$ and
${\bm{u}}^{(i)}({\bm{X}}^{(i)})=\bar{{\bm{u}}}^{(i)}({\bm{X}}^{(i)})$. The
surface traction is prescribed over $\Gamma_{t}^{(i)}$:
${\bm{t}}^{(i)}({\bm{X}}^{(i)})=\bar{{\bm{t}}}^{(i)}({\bm{X}}^{(i)})$. The
surface $\Gamma_{c}^{(i)}$ denotes regions that will potentially contact with
the other body. In the deformed configurations $\omega^{(i)}$ with boundary
$\partial\omega^{(i)}$, these three surfaces occupy regions
$\partial\omega^{(i)}$, $\gamma_{u}^{(i)}$, $\gamma_{t}^{(i)}$ and
$\gamma_{c}^{(i)}$, respectively.

Figure 1: Reference and deformed configurations of bodies in potential
contact.

The unilateral contact constraint is based on the principle of impenetrability
of matter, i.e.,

$\displaystyle\omega^{(1)}\cap\omega^{(2)}=\emptyset.$ (1)

We impose (1) via the gap function $g:\gamma_{c}^{(2)}\to\mathbb{R}$ and
constraint

$\displaystyle g({\bm{x}}^{(2)})\geq 0,\
\forall{\bm{x}}^{(2)}\in\gamma_{c}^{(2)}.$ (2)

To define $g$, we assume that each location
${\bm{x}}^{(2)}\in\gamma_{c}^{(2)}$ has a unique outward unit normal
${\bm{n}}^{(2)}({\bm{x}}^{(2)})$. The projection along
${\bm{n}}^{(2)}({\bm{x}}^{(2)})$ is used to identify a potential contact point
${\bm{x}}^{(1)}\in\gamma_{c}^{(1)}$. The signed distance between these points
define the gap, i.e., $g={\bm{n}}^{(2)}\cdot({\bm{x}}^{(1)}-{\bm{x}}^{(2)})$
[24, 22]. The constraint (2) is one-sided and can bias the analysis [25], but
in general is effective and widely adopted.

To simplify the ensuing presentation, we define
${\bm{u}}=({\bm{u}}^{(1)},{\bm{u}}^{(2)})$,
${\bm{X}}=({\bm{X}}^{(1)},{\bm{X}}^{(2)})$ and
${\bm{x}}=({\bm{x}}^{(1)},{\bm{x}}^{(2)})$. As such, the frictionless contact
for elastic bodies can be posed as an optimization problem:

$\displaystyle\underset{\begin{subarray}{c}{\bm{u}}\end{subarray}}{\text{minimize}}$
$\displaystyle f({\bm{F}},{\bm{u}},{\bm{X}})$ (3) subject to $\displaystyle
g({\bm{x}})\geq 0,\ \forall{\bm{x}}\in\gamma_{c}^{(2)},$

where $f$ is the total energy, i.e.,

$\displaystyle
f({\bm{F}},{\bm{u}},{\bm{X}})=\sum_{i=1}^{2}\int_{\Omega^{(i)}}\phi^{(i)}_{e}({\bm{F}}^{(i)},{\bm{X}}^{(i)})dv+\int_{\Omega^{(i)}}\phi_{ext}^{(i)}({\bm{u}}^{(i)},{\bm{X}}^{(i)}).$
(4)

We apply standard Galerkin finite element discretization to where in the
displacement field is discretized with piecewise smooth interpolation
functions $\Phi^{(i)}_{j}:\Omega^{(i)}\to\Omega$, local coordinates
${\bm{\xi}}\in\mathbb{R}^{3}$, nodal displacement vectors
${\bm{u}}_{h,j}^{(i)}$ and ${\bm{u}}_{h,j}$ such that

$\displaystyle{\bm{u}}({\bm{\xi}},{\bm{u}}_{h,j})=\sum_{j=1}^{n}\Phi_{j}({\bm{\xi}}){\bm{u}}_{h,j},$
(5)

where $n$ is the number of nodes. The subscript $h$ is used to denote nodal
vectors and other discretized quantities, e.g., $\gamma_{h}^{(i)},i=1,2$. The
discretized energy objective for an elastic material is

$\displaystyle
f_{h}({{\bm{u}}}_{h})=\sum_{i=1}^{2}\frac{1}{2}{\bm{u}}_{h}^{(i),T}K^{(i)}{\bm{u}}_{h}^{(i)}-{\bm{f}}_{ext}^{(i),T}{\bm{u}}_{h}^{(i)}=\frac{1}{2}{\bm{u}}_{h}^{T}K{\bm{u}}_{h}-{\bm{f}}_{ext}^{T}{\bm{u}}_{h}.$
(6)

where $K$ is the stiffness matrix and an integrated quantity, and
${\bm{f}}_{ext}$ is the external force vector.

There are multiple ways to discretize and compute the gap function. We choose
the mortar method which integrates the pointwise gap function $g$ on
$\gamma_{c,h}^{(2)}$. Let $S$ be the index set of all nodes on the discretized
contact boundary $\gamma_{c,h}^{(2)}$, the mortar gap function is defined as a
vector function ${\bm{g}}_{h}:\mathbb{R}^{n}\to\mathbb{R}^{m},m=|S|$, with
components

$\displaystyle{\bm{g}}_{h,j}({\bm{u}}_{h})=\int_{\gamma_{c,h}^{(2)}}\Phi_{j}({\bm{\xi}})g({\bm{\xi}},{\bm{u}}_{h})d{\bm{\xi}},$
(7)

for all $j\in S$. Note that the local coordinates reduce to two dimensions on
the surface. Therefore, the discretized optimization problem of (3) with
mortar contact formulation is

$\displaystyle\underset{\begin{subarray}{c}{\bm{u}}_{h}\in\mathbb{R}^{n}\end{subarray}}{\text{minimize}}$
$\displaystyle f_{h}({\bm{u}}_{h})$ (8) subject to
$\displaystyle{\bm{g}}_{h}({\bm{u}}_{h})\geq\mathbf{0}.$

The KKT conditions of (8) is that there exist Lagrange multipliers
${\bm{\lambda}}_{h}\in\mathbb{R}^{m}$ such that

$\displaystyle\nabla
f_{h}({\bm{u}}_{h})-{\bm{\lambda}}_{h}\nabla{\bm{g}}_{h}({\bm{u}}_{h})=0$ (9)
$\displaystyle{\bm{g}}_{h}({\bm{u}}_{h})\geq 0,\ {\bm{\lambda}}_{h}\geq 0,\
\lambda_{h,i}g_{h,i}=0,\ i=1,...,m.$

Equations (9) comprise the standard contact finite element problem [5] and is
the basis for further deriving sensitivities needed for gradient-based design
optimization. The Lagrange multiplier ${\bm{\lambda}}_{h}$ is often
interpreted as the pressure on the contact surface [32].

##  3 Design optimization and sensitivities

In this section, we consider design optimization problems that are subjected
pressure constraints. Let ${\bm{\rho}}\in\mathbb{R}^{p}$ be the design
variable. The mathematical form of our target problems is

$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}\end{subarray}}{\text{minimize}}$
$\displaystyle a({\bm{\rho}},\bar{{\bm{u}}},\bar{{\bm{\lambda}}})$ (10)
subject to
$\displaystyle{\bm{c}}({\bm{\rho}},\bar{{\bm{u}}},\bar{{\bm{\lambda}}})\geq\mathbf{0},$

where
$a({\bm{\rho}},\cdot,\cdot):\mathbb{R}^{p}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$
and
${\bm{c}}({\bm{\rho}},\cdot,\cdot):\mathbb{R}^{p}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{q}$
are the objective and constraint functions, respectively. The displacement
$\bar{{\bm{u}}}$ and contact pressure $\bar{{\bm{\lambda}}}$ are obtained from
the solutions to a contact mechanics problem (3).

Applying the same finite element discretization, the parameterized and
discretized forms of (10) can be written as

$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}\end{subarray}}{\text{minimize}}$
$\displaystyle a_{h}({\bm{\rho}},\bar{{\bm{u}}}_{h},\bar{{\bm{\lambda}}}_{h})$
(11) subject to
$\displaystyle{\bm{c}}_{h}({\bm{\rho}},\bar{{\bm{u}}}_{h},\bar{{\bm{\lambda}}}_{h})\geq\mathbf{0},$

where given ${\bm{\rho}}$,

$\displaystyle\bar{{\bm{u}}}_{h}=$
$\displaystyle\underset{\begin{subarray}{c}{\bm{u}}_{h}\end{subarray}}{\text{argmin}}\
$ $\displaystyle f_{h}({\bm{u}}_{h},{\bm{\rho}})$ (12) subject to
$\displaystyle{\bm{g}}_{h}({\bm{u}}_{h},{\bm{\rho}})\geq\mathbf{0},$

and $\bar{{\bm{\lambda}}}_{h}$ is the Lagrange multiplier associated with
$\bar{{\bm{u}}}_{h}$. In many engineering applications such as shape
optimization, ${\bm{X}}_{h}$ can be a function of ${\bm{\rho}}$, which results
in (12) being dependent on ${\bm{X}}_{h}$ implicitly.

In order to use gradient-based optimization methods to solve (11), we need to
compute the state sensitivities against ${\bm{\rho}}$ given by

$\displaystyle\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}},\
\frac{\partial\bar{{\bm{\lambda}}}_{h}}{\partial{\bm{\rho}}}.$ (13)

For a low-dimensional design space where the dimension $p$ is small, a primal
method suffices, while adjoint methods can be used when $p$ is large. Denote
the derivative of $f_{h}(\cdot,\cdot)$ and ${\bm{g}}_{h}(\cdot,\cdot)$ with
respect to the $i$th variable by $\partial_{i}f_{h}(\cdot,\cdot)$ and
$\partial_{i}{\bm{g}}_{h}(\cdot,\cdot)$, respectively, where $i=1,2$. The
solutions $\bar{{\bm{u}}}_{h},\bar{{\bm{\lambda}}}_{h}$ satisfy the
parameterized KKT conditions

$\displaystyle\partial_{1}f_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})-\bar{{\bm{\lambda}}}_{h}\partial_{1}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})=\mathbf{0},$
(14)
$\displaystyle{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})\geq\mathbf{0},\
\bar{{\bm{\lambda}}}_{h}\geq\mathbf{0},\
\bar{{\bm{\Lambda}}}_{h}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})=\mathbf{0},$

where $\bar{{\bm{\Lambda}}}_{h}$ is a diagonal matrix whose values are the
vector $\bar{{\bm{\lambda}}}_{h}$.

Next, denoting the second-order derivatives of a function by
$\partial^{2}_{i,j},i,j=1,2$ and differentiating the optimality conditions
(14) with ${\bm{\rho}}$, we have

$\displaystyle\partial_{1,1}^{2}f_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}}+\partial_{1,2}^{2}f_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})-[\partial_{1}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})]^{T}\frac{\partial\bar{{\bm{\lambda}}}_{h}}{\partial{\bm{\rho}}}-$
(15)
$\displaystyle\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}}^{T}[\partial_{1,1}^{2}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})]^{T}\bar{{\bm{\lambda}}}_{h}-[\partial_{1,2}^{2}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})]^{T}\bar{{\bm{\lambda}}}_{h}=\mathbf{0}$
,
$\displaystyle\frac{\partial\bar{{\bm{\Lambda}}}_{h}}{\partial{\bm{\rho}}}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})+\bar{{\bm{\Lambda}}}_{h}\partial_{1}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}}+\bar{{\bm{\Lambda}}}_{h}\partial_{2}{\bm{g}}_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})=\mathbf{0}$
.

Solving the linear system of (15) gives us (13). We can then solve (10) with
gradient-based optimization methods by treating $\bar{{\bm{u}}}_{h}$ and
$\bar{{\bm{\lambda}}}_{h}$ as implicit functions of ${\bm{\rho}}$. The
gradients and Jacobian of the objective and constraint functions with respect
to ${\bm{\rho}}$ are computed with

$\displaystyle\frac{da_{h}}{d{\bm{\rho}}}=\frac{\partial
a_{h}}{\partial{\bm{\rho}}}+\frac{\partial
a_{h}}{\partial\bar{{\bm{u}}}_{h}}\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}}+\frac{\partial
a_{h}}{\partial\bar{{\bm{\lambda}}}_{h}}\frac{\partial\bar{{\bm{\lambda}}}_{h}}{\partial{\bm{\rho}}},$
(16)
$\displaystyle\frac{d{\bm{c}}_{h}}{d{\bm{\rho}}}=\frac{\partial{\bm{c}}_{h}}{\partial{\bm{\rho}}}+\frac{\partial{\bm{c}}_{h}}{\partial\bar{{\bm{u}}}_{h}}\frac{\partial\bar{{\bm{u}}}_{h}}{\partial{\bm{\rho}}}+\frac{\partial{\bm{c}}_{h}}{\partial\bar{{\bm{\lambda}}}_{h}}\frac{\partial\bar{{\bm{\lambda}}}_{h}}{\partial{\bm{\rho}}}.$

The second-order derivatives such as
$\partial_{1,2}^{2}f_{h}(\bar{{\bm{u}}}_{h},{\bm{\rho}})$ are arrived at
directly through the chain rule. In many cases, ${\bm{\rho}}$ could affect
${\bm{X}}_{h}$ and the calculation of the second-order derivatives of the gap
function is obtained again through chain rule.

To solve the contact problem (12), we employ the interior-point solver Ipopt
and HiOp [33, 21], also referred to as the forward solver, with the full
analytical Hessian. To solve the design optimization problem (11) using
gradient-based optimization methods, the sensitivities are computed by solving
the linear system (15) and the derivatives in (16) are calculated. For a
robust and efficient outcome, we use Ipopt with its Hessian approximation
option. The gradient-based optimization is implemented in MATLAB.

##  4 Constrained Bayesian optimization

For the gradient-free optimization approach, we design and apply constrained
Bayesian optimization with Gaussian Process (GP). Bayesian optimization is an
iterative black-box optimization method that does not require the analytic
forms of the objective functions or their derivatives. At each iteration, GP
surrogate models replace both the objective and constraint functions and an
acquisition function is used to select new sample points. Applied to design
optimization, each sample point requires a solution to the parameterized
contact problem (12). The optimization variables are typically bounded by a
closed and bounded set ${\bm{\rho}}\in C$.

A Gaussian process surrogate model assumes a multivariate joint Gaussian
distribution between the variables ${\bm{\rho}}$ and an output function
$a_{h}({\bm{\rho}})$, written as

$\displaystyle\begin{bmatrix}a_{h}({\bm{\rho}}_{1})\\\ \vdots\\\
a_{h}({\bm{\rho}}_{T})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m({\bm{\rho}}_{1})\\\
\vdots\\\
m({\bm{\rho}}_{T})\end{bmatrix},\begin{bmatrix}k({\bm{\rho}}_{1},{\bm{\rho}}_{1})\dots
k({\bm{\rho}}_{1},{\bm{\rho}}_{T})\\\ \vdots\\\
k({\bm{\rho}}_{T},{\bm{\rho}}_{1})\dots
k({\bm{\rho}}_{T},{\bm{\rho}}_{T})\end{bmatrix}\ \right).$ (17)

Here, $\mathcal{N}$ is the normal distribution with $T$ samples. The design
variable samples are denoted as ${\bm{\rho}}_{1},\ldots,{\bm{\rho}}_{T}$. The
function $m:\mathbb{R}^{p}\to\mathbb{R}$ is a user-chosen mean function for
the distribution. The function
$k:\mathbb{R}^{p}\times\mathbb{R}^{p}\to\mathbb{R}$ denotes the covariance
function. The posterior Gaussian probability distribution at a new sample
point ${\bm{\rho}}$ can be inferred using Bayes’ rule:

$\displaystyle
a_{h}({\bm{\rho}})|a_{h}({\bm{\rho}}_{1:T})\sim\mathcal{N}(\mu({\bm{\rho}}),\sigma^{2}({\bm{\rho}}))$
(18) $\displaystyle\mu({\bm{\rho}})\ =\
k({\bm{\rho}},{\bm{\rho}}_{1:T})k({\bm{\rho}}_{1:T},{\bm{\rho}}_{1:T})^{-1}\left(a_{h}({\bm{\rho}}_{1:T})-m({\bm{\rho}}_{1:T})\right)+m({\bm{\rho}}_{1:T})$
$\displaystyle\sigma^{2}({\bm{\rho}})\ =\
k({\bm{\rho}},{\bm{\rho}})-k({\bm{\rho}},{\bm{\rho}}_{1:T})k({\bm{\rho}}_{1:T},{\bm{\rho}}_{1:T})^{-1}\sigma({\bm{\rho}}_{1:T},{\bm{\rho}})\
,$

where the vector ${\bm{\rho}}_{1:T}$ is the notation for
$[{\bm{\rho}}_{1},\dots,{\bm{\rho}}_{T}]$ and

$\displaystyle
k({\bm{\rho}}_{1:T},{\bm{\rho}}_{1:T})=[k({\bm{\rho}}_{1}),\dots,k({\bm{\rho}}_{T});\dots;k({\bm{\rho}}_{T},{\bm{\rho}}_{1}),\dots,k({\bm{\rho}}_{T},{\bm{\rho}}_{T})].$
(19)

The functions $\mu$ and $\sigma^{2}$ are the posterior mean and variance,
respectively. It is common to set the mean function $m$ to constant $0$ [8].
The covariance functions, or kernels, of the GP have significant impact on the
accuracy of the surrogate model. A widely adopted kernel, the Squared
Exponential Covariance Function, also known as the power exponential kernel,
is defined as

$\displaystyle k({\bm{\rho}},{\bm{\rho}}^{\prime};\theta)\ =\
\text{exp}\left(-\frac{\left\|{\bm{\rho}}-{\bm{\rho}}^{\prime}\right\|^{2}}{\theta^{2}}\right)\
,$ (20)

where $\theta$ denotes the hyperparameters of the kernel. The hyperparameters
of the GP surrogate model often are optimized during the training or fitting
step by maximizing the log-marginal-likelihood with an optimization algorithm.

The choice of the acquisition function is critical to an efficient and
effective Bayesian optimization algorithm. If the acquisition function
emphasizes on exploration, then typically sample points are chosen to reduce
variance of the model. Otherwise, sample points are chosen to improve
objective function value. One popular acquisition function is the expected
improvement (EI) [1], whose form can be written as

$\displaystyle EI({\bm{\rho}})\ =\
\begin{cases}(\mu({\bm{\rho}})-a_{h}({\bm{\rho}}^{+})-\xi)\Phi(Z)+\sigma(x)\phi(Z),&\text{if}\
\sigma({\bm{\rho}})>0,\\\ 0,&\text{if}\ \sigma({\bm{\rho}})=0,\end{cases}$
(21)

where

$\displaystyle Z\ =\
\begin{cases}\frac{\mu({\bm{\rho}})-a_{h}({\bm{\rho}}^{+})-\xi}{\sigma({\bm{\rho}})},&\text{if}\
\sigma({\bm{\rho}})>0,\\\ 0,&\text{if}\ \sigma({\bm{\rho}})=0\ .\end{cases}$
(22)

In the unconstrained setting, the variable ${\bm{\rho}}^{+}$ in (21) is
defined as

$\displaystyle{\bm{\rho}}^{+}=\underset{\begin{subarray}{c}{\bm{\rho}}_{i}\in{\bm{\rho}}_{1:T}\end{subarray}}{\text{argmax}}\
a_{h}({\bm{\rho}}_{i})$ (23)

of existing $T$ samples, where $C$ is the bounds on ${\bm{\rho}}$. The trade-
off parameter $\xi\geq 0$ affects the trade-off between global search and
local optimization [15]. The functions $\phi$ and $\Phi$ are the probability
density function (PDF) and cumulative distribution function (CDF) of the
standard normal distribution, respectively. The next sample point is chosen by
maximizing the acquisition function. At the $k$th iteration, this leads to

$\displaystyle{\bm{\rho}}_{k}=$
$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}\in
C\end{subarray}}{\text{argmax}}\ EI_{k}({\bm{\rho}}),$ (24)

where $EI_{k}$ is the $EI$ function using samples accumulated up to the $k$th
iteration. In order to solve (24), optimization algorithms including L-BFGS or
random search can be used.

For constrained Bayesian optimization, the choices of acquisition functions
are limited. Fortunately, $EI$ can be extended to constrained expected
improvement $EI_{C}$, by incorporating feasibility penalties to the
improvement function [9]. The ${\bm{\rho}}^{+}$ in (23) is defined as the
feasible optimal solution from the $T$ samples

$\displaystyle{\bm{\rho}}^{+}=$
$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}_{i}\in{\bm{\rho}}_{1:T}\\\
{\bm{c}}_{h}({\bm{\rho}}_{i})\geq 0\end{subarray}}{\text{argmax}}\
a_{h}({\bm{\rho}}_{i}).$ (25)

The expected constrained improvement acquisition function is

$\displaystyle EI_{C}({\bm{\rho}})=PF({\bm{\rho}})EI({\bm{\rho}}),$ (26)

where $PF({\bm{\rho}})=Pr[c_{h}({\bm{\rho}})\geq 0]$ is the probability of the
constraints being satisfied. For multiple inequality constraints, each
constraint function is modeled as an independent GP surrogate model and is
often simplified to be conditionally independent given ${\bm{\rho}}$. That is,
the probability of a feasible ${\bm{\rho}}$ for constraints
${\bm{c}}_{h}({\bm{\rho}})\geq 0$ satisfies

$\displaystyle Pr({\bm{c}}_{h}({\bm{\rho}})\geq
0)=\prod_{i=1}^{q}Pr({\bm{c}}_{h,i}({\bm{\rho}})\geq 0).$ (27)

The constrained Bayesian optimization algorithm is given in Algorithm 1.

Algorithm 1 Constrained Bayesian optimization with GP

1: Choose initial sampling data ${\bm{\rho}}_{0}$ in the bounds.

2: Train GP surrogate models on the initial data.

3: for $i=0,1,2,\dots$ do

4: Evaluate acquisition function and find
${\bm{\rho}}_{i}=\text{argmax}_{{\bm{\rho}}}EI_{C,k}(x)$ within the bounds.

5: Run finite element simulation with design variables ${\bm{\rho}}_{i}$.

6: Compute the objective $a_{h}({\bm{\rho}}_{i})$ and constraints
${\bm{c}}_{h}({\bm{\rho}}_{i})$ based on the simulation result.

7: Retrain the GP surrogate models with the addition of the new sample
${\bm{\rho}}_{i}$, $a_{h}({\bm{\rho}}_{i})$ and
${\bm{c}}_{h}({\bm{\rho}}_{i})$.

8: Solve (11) as needed with the updated surrogate model. Evaluate the
stopping criteria. Exit if satisfied.

9: end for

The training of the surrogate model refers to finding the optimal
hyperparameters given the sample points. Line $8$ of Algorithm 1 can either
return the feasible sample with the smallest objective or with the smallest
posterior mean. Line $4$ can be solved by interior-point L-BFGS methods, our
default choice, or random search algorithms. One stopping criteria is that the
algorithm exists upon reaching a prescribed maximum number of iterations.

Algorithm 1 is implemented in MATLAB. We note that Bayesian optimization
algorithms are often coded as maximization algorithms. In such case, $-a_{h}$
is the maximization objective.

##  5 Numerical examples

To demonstrate and compare the effectiveness of the gradient-based and
gradient-free optimization methods, we present two numerical examples that are
inspired by real-life engineering problems. These model problems can be easily
extended to more complicated and larger-scale problems using the same
optimization method we presented in this paper.

###  5.1 Wedge problem

The first numerical example is derived from the problem of designing high-
current joints, whose illustration is given in Figure 2. In this case, we
consider a high-current joint where the current and magnetic field diffuse
from the outside (positive z) inwards. At very early times, the current is
concentrated at the very outer fiber of the joint, so it is important large
contact pressure be maintained in the very beginning of the contact region,
i.e., $\bar{{\bm{\lambda}}}_{h}>{\bm{\lambda}}_{l}$, so that arcing and/or
joint separation dot occur. As the current builds up, it diffuses into the
surface to a depth $s=\alpha D$, but it is still important that the contact
pressure be maintained above a minimum value ${\bm{\lambda}}_{l}$ at the outer
fiber and all the way down to the depth $s=\alpha D$, since every additional
increment in current diffuses in from the outside. At the same time, as
current builds, a magnetic pressure is created, represented by $P_{2}(t)$ in
Figure 2, which pushes inwards from the outside. The general design goal is to
create a joint geometry which maintains adequate contact pressure within the
joint, but without exceeding any material limits (e.g. locally
$\bar{{\bm{\lambda}}}_{h}<{\bm{\lambda}}_{m}$). At the optimum design, the
amount of stabilized preload, represented by $P_{1}$ in Figure 3, should be a
minimum, as large values of contact preload require additional structure
surrounding the joint. The joint is modeled in plane strain, but the extension
to 3-D or axisymmetric analysis is trivial. Based on the description above,
the reference configuration in the $y-z$ plane can be seen in Figure 2.

Figure 2: Wedge problem in the $y-z$ plane. $P_{1}(t)$ and $P_{2}(t)$ are two
time-dependent loads.

The two external loads are depicted in Figure 3 below.

Figure 3: The pressure loads $P_{1}(t)$ and $P_{2}(t)$ as they change over
time. The scalars $P_{1}$ and $P_{2}$ are the stabilized values of the loads.

The two wedges are subjected to two external loads in the form of pressure,
where the first pressure $P_{1}(t)$ ramps up linearly from time $0$ and value
$0$ to $5$ seconds and value $P_{1}$. The second pressure $P_{2}(t)$ ramps up
linearly from $5$ to $10$ seconds for a value from $0$ to $2$. The wedges are
restricted in their movement in the $x$-direction and the left wedge has zero
$z$-displacement on its left-most surface. The material for both wedges are
chosen to be linear isotropic elastic material with the elastic modulus
$E=200$ and Poisson’s ratio $\nu=0.3$. For demonstration purpose, the problem
is run unitless and can be easily extended to more complicated materials. Both
wedges are discretized with linear hexahedral elements, as shown in Figure 4.

Figure 4: The reference configurations of the wedges in the initial design.

Given the plane strain boundary conditions imposed, the number of elements in
the $x$ direction is reduced to $1$ to improve computational efficiency. The
contact surface on the right wedge is chosen to be the surface for the mortar
contact integrals.

As mentioned earlier, the design objective the wedge problem is the stabilized
pressure $P_{1}$, while we impose constraints on the contact pressure on the
upper part of the contact surfaces at both $t=5$ and $t=10$ seconds. The
mathematical expression of the problem is

$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}\in\mathbb{R}^{3}\end{subarray}}{\text{minimize}}$
$\displaystyle P_{1}$ (28) subject to
$\displaystyle\lambda_{i}^{e}\geq\lambda_{l},\ s\leq\alpha D,$
$\displaystyle\lambda_{m}\geq\lambda_{i}^{e},\ \forall s,$

where the parameters are set to $\lambda_{m}=20$, $\lambda_{l}=1.0$ and
$\alpha=0.4$ so that the top part of the contact surface are required to be in
active contact. To reduce the potential nonsmoothness of pressure, the
constraints are imposed on element segment pressure $\lambda_{i}^{e}$ for
segment $i$, computed through the average of nodal pressure. For the given
discretization, we choose an element segment made of two elements, with a
total of $11$ element segments, $22$ elements and $23$ nodes on the contact
surface of the right wedge. Therefore, the constraints in (28) are imposed on
a total of four element segments, equivalent to eight elements, towards the
top of the right wedge. The optimization variables are the two angles of the
wedges $\theta_{1}$, $\theta_{2}$ and the pressure $P_{1}$. All of the three
variables have bound constraints to ensure contact would occur at the two time
steps. Specifically, the bound constraints are

$\displaystyle 30^{\circ}\leq\theta_{1},\theta_{2}\leq 60^{\circ},$ (29)
$\displaystyle 0.5\leq P_{1}\leq 1.5.$

The initial values of the variables are set to $\theta_{1}=39^{\circ}$,
$\theta_{2}=41^{\circ}$ and $P_{1}=1$. Using the gradient-based approach, the
optimal solution is successfully found in 35 iterations with the values
$\theta_{1}=30,\theta_{2}=31.9583$ and $P_{1}=0.7527$. The deformed
configuration of the optimal design at time $5$s and $10$s are shown in Figure
5.

Figure 5: The deformed configurations of the optimal design at time 5$s$ (top)
and 10$s$ (bottom).

The optimization paths are plotted in Figure 6.

Figure 6: Optimization result for the wedge problem. The objective and
constraint violation values are plotted on the left. The logarithm dual
optimality measure is given on the right.

Next, we report the results from constrained Bayesian optimization described
in section 4. We use Latin hypercube sampling to generate $8$ initial samples
for the three dimensional design space. One more feasible sample point is
added to start the algorithm. The algorithm is given a total budget of $50$
iterations, comparable to the number of iterations required for convergence in
the gradient-based approach. Thanks to the randomness of the algorithm, we
repeat the same run $10$ times and present the average results.

Since the feasible region is quite small, as evidenced by the high ratio of
infeasible iterates during the optimization in Figure 6, many samples are
infeasible. However, the constrained Bayesian optimization succeeded in
finding improved solution in $50$ iterations compared to the initial design.
Figure 7 shows the objective and constraint violation of the samples collected
during one of the runs, where the feasible and infeasible ones are marked with
different signs. The optimal design variables are given on the right in Figure
7, where the mean value and standard deviation of the ten constrained Bayesian
optimization runs are computed and plotted, together with the optimal solution
obtained through the gradient method.

Figure 7: Constrained Bayesian optimization result for the wedge problem.

We note that throughout all the iterations and runs, the design obtained
through gradient-based approach remains the optimal one.

The optimization lays heavy emphasis on contact pressure, which is a critical
source of nonsmoothness as discussed in previous sections. Hence, to better
illustrate the feasibility region and how the pressure depends on optimization
variables, we plot the pressure distribution of the $23$ nodes on the contact
surface of the right wedge, as well as the averaged contact pressure of the
$11$ element segments. The plot for the initial design
$\theta_{1}=39,\theta_{2}=41,P_{1}=1$ is given in Figure 8. The $x$-axis
represents the nodes on the contact surface on the $x=0$ plane. The order of
increasing index matches with that of a decreasing $y$ coordinate, i.e., node
$1$ is the upper-most node and element segment $1$ is the upper-most element
segment on the contact surface. The upper and lower bounds of pressure on the
first four element segments are marked on the plots as well and the $y$-axis
uses symlog scale for better illustration.

Figure 8: Pressure distribution on the contact surface for the initial design:
nodal pressure (left) and averaged pressure on element segments (right).

The pressure plot for the optimal design
$\theta_{1}=30,\theta_{2}=31.9583,P_{1}=0.7527$ is given in Figure 9. It is
clear that the initial design does not produce a desired pressure
distribution. Moreover, the averaging of pressure over element segments
significantly smooths the nodal pressure distribution, which possibly leads to
better convergence behavior of gradient-based optimization algorithms.
Finally, we note that the two time steps of concern typically have quite
different pressure distribution thanks to the different loading conditions.
Therefore, satisfying pressure constraints at both time steps prove to be
challenging.

Figure 9: Pressure distribution on the contact surface for the optimal design:
nodal pressure (left) and averaged pressure on element segments (right).

###  5.2 Marman clamp problem

The second example considers the design of a Marman clamp [3]. These devices
are commonly used to connect two cylindrical sections together, see Figures 10
and 11. In a general Marman clamp, the cylinders terminate in flanges which
meet at a plane. The two flanges interface with a retainer, which is forced
against the flanges by the action of a clamp. The retainers are commonly
manufactured in segments to allow for inwards displacement without generating
significant hoop stresses. The shape of the retainer and the flanges forces
the flanges together. In [34], the authors optimized clamp band joint using
the so-called constrained sequential approximate optimization method, which
bears some similarity to the constrained Bayesian optimization method. In [4],
uncertainty analysis for the clamp band joints is incorporated into a
mathematical model.

In the manifestation considered here, the clamp is in the form of a band which
is forced against the retainer by designing the components such that there is
an interference fit between them. The hoop stress generated in the band is
reacted by a contact pressure between the band and the retainer, which is
further reacted by contact pressure between the retainer and the flanges, and
ultimately the two flanges are forced together, inducing contact pressure
between them.

Figure 10: View looking down the cylindrical axis of a generic Marman clamp.
The flange is part of the cylindrical object that is being clamped together.
The gap between the segments of the retainer ensures that hoop stresses in the
retainer are minimized. The band is either stretched over the retainer such
that the induced hoop stresses in the band induce compressive force on the
retainer, or the compression may be induced by splitting the band at some
point and pulling the pieces together using a screw mechanism.

The cylindrical structures can be treated as axisymmetric so that only a small
portion of the clamp is simulated and the angular coordinates ($\theta$) are
fixed in the cylindrical coordinate system. Given a $\theta$, the bodies in
the radial and longitudinal plane, $x-y$ plane, is shown in Figure 11.

Figure 11: Left: Cross section of the clamp through a plane of constant
$\theta$. The cylinder is split into two pieces, the walls of the cylinder are
terminated by flanges which meet at the centerline of the clamp. The retainer
is forced against the flanges by the band, causing the flanges to be clamped
against each other. Right: the contact surfaces and boundary conditions on an
example design. The pressure involved in the constraints are marked.

Similarly, only the top flange (referring to Figure 11) and the top half of
the band and the clamp are simulated. Three contact planes are defined –
contact between the band and the retainer, between the retainer and the top
flange, and between the top flange and the vertical symmetry plane,
represented by an artificial box which is held rigid. The vertical
displacement at the vertical symmetry plane of the retainer and the band is
held fixed in accordance with typical symmetry boundary conditions. As
mentioned earlier, hoop stress in the band is generated through an initial
interference fit, represented as an initial overlap in the meshes.

Both the flange and the band are given isotropic linear elastic materials with
the elastic modulus at $E=2\times 10^{6}$ and Poisson’s ratio at $\nu=0.3$.
The retainer is given an anisotropic material that has a smaller ($0.1\%$)
elastic and shear modulus in the $\theta$-direction in order to simulate the
gap between segments of the retainer (Figure 10).

The referential three-dimensional finite element model of an example design
(the initial design) is shown on the left in Figure 12. The mesh around the
potential contact surface between the retainer and the flange is refined so
that a more accurate contact solution can be obtained. Notice that the initial
configurations have penetrations, which are resolved as part of the contact
forward solution.

Figure 12: Finite element mesh of the initial design in the undeformed
configuration (left) and deformed configuration (right). Note that the band is
meshed with an initial penetration which will be forced out by the contact,
inducing hoop stresses in the band and compressive stresses in the retainer.
The bodies are meshed as a 3D wedge with angled boundary conditions to enforce
axisymmetry.

To ensure robustness of the design against separation of the flange, the
objective is chosen to be minimization of the compliance, that is maximizing
the stiffness of the overall system. Meanwhile, the maximum contact pressure
between the retainer and the flange is constrained by user-defined upper and
lower bounds. The maximum pressure is imposed on the element pressure to
reduce the nonsmoothness of the nodal pressure distribution. The element
pressure is denoted as $\bar{{\bm{\lambda}}}^{e}$ and computed as the average
of its nodal pressure . The $40$ elements on the contact surface are used to
define the constraints. To prevent leaks from the cylindrical assembly and
ensure adequate sealing, we enforce an interface force at the bottom-left part
of the flange. This is reflected as a constraint on the sum of the four nodal
pressure at the bottom-left part of the flange, given that the mesh and
surface area of this region is independent of the design variable.

The design variables are the contact geometries of the retainer and flange.
Specifically, in the $x-y$ plane the contact surfaces are defined by two cubic
Bézier curves, on both the retainer and the flange, an example of which is
shown in Figure 13. Given four nodes ${\bm{P}}_{i},i=1,2,3,4$, the explicit
form of a curve is

$\displaystyle{\bm{B}}(t)=(1-t)^{3}{\bm{P}}_{1}+3(1-t)^{2}t{\bm{P}}_{2}+3(1-t)t^{2}{\bm{P}}_{3}+t^{3}{\bm{P}}_{4},\
0\leq t\leq 1,$ (30)

where the vector ${\bm{B}}\in\mathbb{R}^{3}$ is the nodal coordinate. The
curve passes through ${\bm{P}}_{1}$ and ${\bm{P}}_{4}$ but not necessarily the
two middle control points ${\bm{P}}_{2}$ and ${\bm{P}}_{3}$.

Given the geometry and assembly constraints, for both the retainer and the
flange, the two end points ${\bm{P}}_{1}$ and ${\bm{P}}_{4}$ of the two curves
are fixed, as well as the $x$-coordinates of the two middle control points
${\bm{P}}_{2}$ and ${\bm{P}}_{3}$ of each curve. The $y$-coordinates of the
four middle control points, two from each curve, are the optimization
variables. Their $z$-coordinates are determined by cylindrical symmetry, where
$\theta=2^{\circ}$ is fixed. The inner radius of the cylinder as seen in
Figure 10 is $20$.

Figure 13: The curve to be optimized in the $x-y$ plane. The $y$-coordinates
of the ‘variable’ nodes are design variables.

Further, each variable is bounded in order to ensure valid retainer and flange
shapes as well as contact between the retainer and flange. The bounds are
simply chosen through numerical experiments and observation.

The mathematical form of the optimization problem is

$\displaystyle\underset{\begin{subarray}{c}{\bm{\rho}}\in\mathbb{R}^{4}\end{subarray}}{\text{minimize}}$
$\displaystyle a({\bm{\rho}},\bar{{\bm{u}}}_{h})$ (31) subject to
$\displaystyle\lambda_{l,o}\leq\bar{\lambda}_{h,o}$
$\displaystyle\bar{\lambda}_{l}\leq\max(\bar{{\bm{\lambda}}}^{e})\leq\bar{\lambda}_{u}$
$\displaystyle{\bm{\rho}}_{l}\leq{\bm{\rho}}\leq{\bm{\rho}}_{u},$

where ${\bm{\rho}}\in\mathbb{R}^{4}$, $\bar{{\bm{\lambda}}}^{e}$ is the
pressure vector of the element segments on the contact surface between the
flange and the retainer. The objective function $a$ is the compliance
function. The user-defined constraint bounds heavily influence the optimal
design and can be tuned for other design objectives. The pressure
$\bar{\lambda}_{h,o}$ is the sum of pressure of the four nodes on the bottom-
left part of the flange. In this example, we set the $x$-coordinates of the
four middle control points to $[0.617,0.883]$ for the retainer and
$[0.50,0.773]$ for the flange so that all four control points are evenly
distributed in the $x$-direction. We choose the bounds to be
${\bm{\rho}}_{u}=[0.408,0.354,0.47,0.44]$ and
${\bm{\rho}}_{l}=[0.35,0.35,0.443,0.3834]$ and constraint parameters to be
$\lambda_{l,o}=30$, $\bar{\lambda}_{l}=300$ and $\bar{\lambda}_{u}=650$.

The initial ${\bm{\rho}}$ is set to $[0.35,0.32,0.45,0.41]$. Its deformed
configuration from contact simulation is shown on the right in Figure 12. The
pressure on the bottom-left nodes of the flange is essentially $0$, as the
contact between the flange and the retainer occurs towards the right end of
the flange. For the gradient-based approach, an optimal solution is
successfully found in 48 iterations with Ipopt, using the adaptive barrier
update strategy to a dual optimality accuracy of $10^{-4}$. The optimal set of
values are $[0.40713,0.35344,0.443,0.38343]$, with three of the four variables
not reaching bounds.

The objective, constraint violation and dual optimality during the
optimization are plotted in Figure 14, where the $y$-axis of constraint
violation uses symlog scale for better illustration. Additionally, the
feasibility restoration steps are marked with the $\times$ sign.

Figure 14: The optimization iterations for the clamp problem using Ipopt. The
objective and constraint violation values are given on the left. The logarithm
of dual optimality measure is given on the right.

The deformed configuration of the optimal design is shown in Figure 15.
Compared to the initial design, the contact region has moved to the left so
that pressure could be passed onto the bottom-left part of the flange.

Figure 15: The undeformed (left) and deformed (right) configurations of the
optimal design.

The pressure distribution of the bottom contact surface and the contact
surface between the flange and the retainer is plotted in Figure 16. For the
contact surface between the flange and the retainer, we plot the element
pressure. It is worth noting that thanks to the curved design of the two
contacting surfaces, the element pressure on the contact surface between the
flange and the retainer is not smooth or monotonic.

Figure 16: Pressure distribution of the initial and optimal design at the
bottom contact surface (left) and the contact surface between the flange and
the retainer. The node index for the bottom contact surface increases with the
$x$-coordinate. The element index of the contact surface between the flange
and the retainer increases with the $x$-coordinate as well.

While the contact nodes and regions change with the design variables,
sensitivities with Ipopt manage to overcome potential nonsmoothness to find a
reasonable solution through careful design of the bound constraints. Future
work to better understand the nonsmoothness and design more robust
optimization algorithms is being planned.

Next, we apply the constrained Bayesian optimization algorithm (Algorithm 1)
to the clamp design problem, where $16$ initial Latin hypercube samples are
used due to the increased dimension of the design space. One more feasible
sample point is added to start the algorithm. We run the Bayesian algorithm
for $50$ iterations and repeat the same run $10$ times to obtain less
uncertain results.

As evidenced by the large number of feasibility restoration steps with Ipopt,
the feasible region for this problem is rather small. Applying BFGS interior-
point method to find the maximum of the constrained expected improvement
acquisition function (line 4 of Algorithm 1) is not effective due to the low
probability of constraints being satisfied. Therefore, we use random search
method for this problem where the maximum constrained expected improvement
value is chosen from $10000$ random variables at each iteration to determine
the next sample. We note that only the constrained expected improvement
function, which has a closed form, is evaluated at these variables. No
additional contact forward simulation is required.

On average, we obtain $6$ out of the $66$ total samples as feasible designs.
We compute the relative absolute error of the mean value of the optimal design
from Bayesian optimization from the optimal design from Ipopt, as plotted in
17. Overall, the Bayesian optimization algorithm generates designs close to
the one from the gradient-based method.

Figure 17: The relative absolute error of the mean optimal variables from
constrained Bayesian optimization. Plot is in log scale.

The overall results from the constrained Bayesian optimization show that an
improved design can be obtained for a reasonable sample size. Compared to the
gradient-based approach, it could offer an alternative way for optimization
with higher efficiency but reduced accuracy, particularly for design problems
with a relatively small number of variables and no available derivatives.
However, the constrained Bayesian optimization algorithm adopted here requires
a feasible initial sample, which could limit its applicability. The authors
are proposing a new and improved constrained Bayesian optimization method in
[31]. We also point out that for the general topology optimization problems
where the densities are the variables, the large dimension of the problem
could make Bayesian optimization impractical.

##  6 Conclusion

In this paper, we solved design optimization problems in unilateral contact
under pressure constraints, using both gradient-based optimization approach
and a gradient-free constrained Bayesian optimization algorithm. The contact
finite element simulation uses the mortar gap function formulation and
interior-point forward solver. Through smoothing of nodal pressure and careful
design of the bound constraints, we successfully find optimal solutions using
the gradient-based method and sensitivities, even though the problems are
considered nonsmooth. The Bayesian optimization approach produced relatively
good results and could be suitable for problems without reliable sensitivity
information. To solve our target problems more efficiently and robustly, we
should address the nonsmoothness directly and develop algorithms that might
overcome it. We aim to work on this challenge in the future.

## Acknowledgments

The authors would like to thank Professor Daniel A. Tortorelli for his help
and support in preparation of this manuscript. The valuable feedback from
Professor Tortorelli has contributed to the manuscript. Prepared by LLNL under
Contract DE-AC52-07NA27344. Release number: LLNL-JRNL-863659-DRAFT.

## References

  * [1] E. Brochu, V. M. Cora, and N. De Freitas.  A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  arXiv preprint arXiv:1012.2599, December 2010.
  * [2] M. Bruggi and P. Duysinx.  A stress–based approach to the optimal design of structures with unilateral behavior of material or supports.  Structural and Multidisciplinary Optimization, 48:311–326, 2013\.
  * [3] NASA Goddard Space Flight Center.  Marman clamp system design guidelines, guideline no. gd-ed2214.  Technical report, 2000.
  * [4] A. Davar, R. Azarafza, and S. Mousavi.  Uncertainty analysis and reliability estimation of clamp band joints.  Materials Science and Technology, 36(13):1487–1502, 2020.
  * [5] C. Eck, J. Jarusek, and M. Krbec.  Unilateral contact problems: variational methods and existence theorems, volume 270.  CRC Press, 2005.
  * [6] E. Fancello.  Topology optimization for minimum mass design considering local failure constraints and contact boundary conditions.  Structural and Multidisciplinary Optimization, 32:229–240, 09 2006\.
  * [7] F. Fernandez, M. A. Puso, J. Solberg, and D. A. Tortorelli.  Topology optimization of multiple deformable bodies in contact with large deformations.  Computer Methods in Applied Mechanics and Engineering, 371:113288, 2020.
  * [8] P. I. Frazier.  Bayesian optimization.  In Recent advances in optimization and modeling of contemporary problems, pages 255–278, October 2018.
  * [9] J. R. Gardner, M. J. Kusner, Z. Xu, K. Q. Weinberger, and J. P. Cunningham.  Bayesian optimization with inequality constraints.  In Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32, ICML’14, pages II–937–II–945. JMLR.org, 2014.
  * [10] J. O. Hallquist, G. L. Goudreau, and D. J. Benson.  Sliding interfaces with contact-impact in large-scale lagrangian computations.  Computer Methods in Applied Mechanics and Engineering, 51(1):107–137, 1985.
  * [11] D. Hilding, A. Klarbring, and J. Petersson.  Optimization of Structures in Unilateral Contact.  Applied Mechanics Reviews, 52(4):139–160, 04 1999.
  * [12] Y. C. Hsu, Y. W. Gung, S. L. Shih, C. K. Feng, S. H. Wei, C. H. Yu, and C. S. Chen.  Using an optimization approach to design an insole for lowering plantar fascia stress—a finite element study.  Annals of biomedical engineering, 36:1345–1352, 2008.
  * [13] H. Kristiansen, K. Poulios, and N. Aage.  Topology optimization for compliance and contact pressure distribution in structural problems with friction.  Computer Methods in Applied Mechanics and Engineering, 364:112915, 2020.
  * [14] T. A. Laursen.  Formulation and treatment of frictional contact problems using finite elements.  Stanford University, 1992.
  * [15] D. J. Lizotte.  Practical Bayesian optimization.  PhD thesis, University of Alberta, Edmonton, Alberta, Canada, 2008.
  * [16] A. Mathern, O. S. Steinholtz, A. Sjöberg, et al.  Multi-objective constrained Bayesian optimization for structural design.  Structural and Multidisciplinary Optimization, 63:689–701, February 2021.
  * [17] C. Niu, W. Zhang, and T. Gao.  Topology optimization of continuum structures for the uniformity of contact pressures.  Structural and Multidisciplinary Optimization, 60:185–210, 2019\.
  * [18] J. Nocedal and S. J. Wright.  Numerical Optimization.  Springer, New York, 2nd edition, 2006.
  * [19] C. B. Pedersen.  Crashworthiness design of transient frame structures using topology optimization.  Computer methods in applied mechanics and engineering, 193(6-8):653–678, 2004.
  * [20] J. Petersson and M. Patriksson.  Topology optimization of sheets in contact by a subgradient method.  International journal for numerical methods in engineering, 40(7):1295–1321, 1997.
  * [21] C. G. Petra, N. Y. Chiang, and J. Wang.  Hiop–user guide.  Center for Applied Scientific Computing, Lawrence Livermore National Labo-ratory, Tech. Rep. LLNL-SM-743591, 2018.
  * [22] M. A. Puso and T. A. Laursen.  A mortar segment-to-segment contact method for large deformation solid mechanics.  Computer Methods in Applied Mechanics and Engineering, 193(6):601–629, 2004.
  * [23] M. A. Puso and T. A. Laursen.  A mortar segment-to-segment frictional contact method for large deformations.  Computer Methods in Applied Mechanics and Engineering, 193(45):4891–4913, 2004.
  * [24] F. Sewerin and P. Papadopoulos.  On the finite element solution of frictionless contact problems using an exact penalty approach.  Computer Methods in Applied Mechanics and Engineering, 368:113108, 2020.
  * [25] J. M. Solberg, R. E. Jones, and P. Papadopoulos.  A family of simple two-pass dual formulations for the finite element solution of contact problems.  Computer Methods in Applied Mechanics and Engineering, 196(4):782–802, 2007.
  * [26] N. Strömberg and A. Klarbring.  Topology optimization of structures in unilateral contact.  Structural and Multidisciplinary Optimization, 41:57–64, 2010.
  * [27] K. Svanberg.  A globally convergent version of mma without linesearch.  In Proceedings of the first world congress of structural and multidisciplinary optimization, volume 28, pages 9–16. Goslar, Germany, 1995\.
  * [28] İ. Temizer, M.M. Abdalla, and Z. Gürdal.  An interior point method for isogeometric contact.  Computer Methods in Applied Mechanics and Engineering, 276:589–611, 2014.
  * [29] J. Wang, N. Y. Chiang, A. Gillette, and J. L. Peterson.  A multifidelity Bayesian optimization method for inertial confinement fusion design.  Physics of Plasmas, 31(3), 2024.
  * [30] J. Wang and P. Papadopoulos.  Optimization of process parameters in additive manufacturing based on the finite element method, 2023.
  * [31] J. Wang, C. G. Petra, and J Luc Peterson.  Constrained Bayesian optimization with merit functions.  arXiv preprint arXiv:2403.13140, 2024.
  * [32] P. Wriggers and T. A. Laursen.  Computational contact mechanics, volume 2.  Springer, 2006.
  * [33] A. Wächter and L. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Math. Program., 106(1):25–27, 2006.
  * [34] B. Yu, Y. Lei, Z. Wang, D. Zhang, and H. Li.  Lightweight design of a clamp band joint via a constrained sequential approximate optimization method in the irregular design domain.  Engineering Optimization, 56(2):155–178, 2024.
  * [35] G. Zavarise and L. De Lorenzis.  The node-to-segment algorithm for 2d frictionless contact: classical formulation and special cases.  Computer Methods in Applied Mechanics and Engineering, 198(41-44):3428–3451, 2009.
  * [36] C. Zillober.  A globally convergent version of the method of moving asymptotes.  Structural optimization, 6:166–174, 1993.