Topology Optimisation under Uncertainties with Neural Networks

Abstract

Topology optimisation is a mathematical approach relevant to different engineering problems where the distribution of material in a defined domain is distributed in some optimal way, subject to a predefined cost function representing desired (e.g., mechanical) properties and constraints. The computation of such an optimal distribution depends on the numerical solution of some physical model (in our case linear elasticity) and robustness is achieved by introducing uncertainties into the model data, namely the forces acting on the structure and variations of the material stiffness, rendering the task high-dimensional and computationally expensive. To alleviate this computational burden, we develop two neural network architectures (NN) that are capable of predicting the gradient step of the optimisation procedure. Since state-of-the-art methods use adaptive mesh refinement, the neural networks are designed to use a sufficiently fine reference mesh such that only one training phase of the neural network suffices. As a first architecture, a convolutional neural network is adapted to the task. To include sequential information of the optimisation process, a recurrent neural network is constructed as a second architecture. A common 2D bridge benchmark is used to illustrate the performance of the proposed architectures. It is observed that the NN prediction of the gradient step clearly outperforms the classical optimisation method, in particular since larger iteration steps become viable.

1. Introduction

Structural topology optimisation is the (engineering-oriented) process of designing a construction part using optimisation algorithms under certain constraints. The resulting designs usually have a large influence on the subsequent production costs. The starting point of the process is a design domain that represents the maximum space available for the optimised component to be developed. The outcome presents information about which parts of the design space are occupied by material and which are void. Often, the task is motivated by mechanical requirements, e.g., sufficient stiffness of the constructed part with respect to assumed forces acting on it while certain predetermined points or surfaces should connect to other parts. A typical physical model for this comes from linear elasticity, describing the displacement field given material properties and forces. For the mathematical optimisation, it is repeatedly necessary to compute the stress distribution determined by the physical model in the design domain (more precisely, in the parts of the domain with material). This potentially complex computational task usually relies on the finite element method (FEM), which is based on a discretisation of the domain into elements. Most commonly, the domain is represented as mesh consisting of disjoint simplices, i.e., triangles in 2D and tetrahedra in 3D.

Since the optimisation process easily requires several hundred evaluations of the state equation to evolve the material distribution, it is of significant interest to develop techniques that reduce this computational burden. This becomes much more pronounced when uncertainties of the model data should be considered in the computations. The treatment of uncertainties has been developed extensively from a theoretical and practical point of view over the last decade in the area of Uncertainty Quantification (UQ). A common way to describe uncertainties is by means of random fields, whose (Karhunen-Loève) expansions depend on a possibly very large number of random variables. The parameter space spanned by these random variables leads to very high-dimensional state problems for which the derived optimisation problems are very difficult to solve.

This paper investigates the application of a trend in scientific computing for current topology optimisation methods, namely the use of modern machine learning techniques. More precisely, our objective is to improve the efficiency of the structural topology optimisation problem by predicting gradient steps based on generated training data. This efficiency gain directly transfers to our ability to compute much more involved risk-averse stochastic topology optimisation problems with random data. In this case, the topology is optimised with an adjusted cost functional including the CVaR (conditional value at risk), by which unlikely events can be taken into account in contrast to simply optimising with the mean value of possible load scenarios. In addition to random loads, we also include random material properties which, e.g., can enter the model in terms of material errors or impurities. We emphasise that risk-averse optimisation based on some risk measure is a timely topic, which plays a role in many application areas. Despite its relevance, this type of problem has not been covered widely in the literature yet. In fact, the authors are not aware of any other machine-learning-assisted approach for risk-averse topology optimisation. This might be due to the more involved mathematical framework and substantially higher computational complexity. To achieve performance benefits with topology optimisation in this paper, we adapt concepts from the field of deep learning to approximate multiple iterations of the optimisation process and render the overall optimisation more efficient.

The goal of topology optimisation is to satisfy the technical requirements of a component (for instance, stiffness with respect to certain loading scenarios) with minimal use of material. There are different approaches to describe the topology in a flexible way such that substantial changes are possible. We follow our previous work in [1] and use a phase field model which describes the density of material with a value in $[0,1]$. The starting point is the definition of a physical design space available for the component under consideration. This space is completely filled with a material in the sense of an initial solution. Furthermore, all points at which loads act on the component, as well as the type of the respective load, are prescribed. The optimisation aims to achieve a homogeneous, minimum possible deformation at all optimised points of the component under the imposed (possibly continuous and thus infinitely many) loading scenarios. Here, a minimum compliance corresponds to a maximum stiffness. In general, even solving the underlying partial differential equation (PDE) of this problem for deterministic coefficients of the PDE presents a complex task. Furthermore, PDE coefficients which describe material and the loads have a strong influence on the resulting topology, i.e., even small changes in these coefficients can lead to large differences in the resulting topology. This results in considerable computational effort in the stochastic settings, since the solution has to be calculated for a sufficient number of data realisations to become reliable. Hence, the modelling of these stochastic settings for example (with the most obvious approach) by a Monte Carlo (MC) simulation increases the required iteration steps linearly in the number of simulations.

A method to numerically tackle topology optimisation uncertainty was presented in [2]. In this paper, we extend the previous work by introducing Deep Neuronal Networks (DNN) that are designed to provide a prediction of the next gradient step. Since topologies discretised with finite elements can be represented as images, Convolutional Neural Networks (CNN) seem natural candidate architectures for this task and there has already been some research on this approach for the deterministic setting. An introduction is presented in [3] where the conventional topology optimisation algorithms are replicated in a computationally inexpensive way. Furthermore, a CNN is used in [4] to approximate the last iteration steps of a gradient method of a topology optimisation after a fixed number of steps to refine a “fuzzy” solution. A CNN architecture is also used in [5] to solve a topology optimisation problem and trained with large amounts of data. The resulting NNs were able to solve problems with boundary conditions different to their training data. In [6], the problem is stated as an image segmentation task and a deep NN with encoder–decoder architecture is leveraged for pixel-wise labeling of the predicted topology. Another encoding–decoding U-Net CNN architecture is presented in [7], providing up- and down-sampling operators based on training with large datasets. In [8] a multilevel topology optimisation is considered where the macroscale elastic structure is optimised subject to spatially varying microscale metamaterials. Instead of density, the parameters of the micromaterial are optimised in the iteration, using a single-layer feedforward Gaussian basis function network as surrogate model for the elastic response of the microscale material.

A discussion on solving PDEs with the help of Neural Networks (NN) for instance of the Poisson equation and the steady Navier–Stokes equations is provided in [9]. In a relatively new approach, through a combination of Deep Learning and conventional topology optimisation, the Solid Isotropic Material with Penalisation (SIMP) was presented in [10], which could reduce the computational time compared to the classical approach. The authors use a similar method as [4] except that the underlying optimisation algorithm performs a mesh refinement after a fixed number of iterations. To improve this step, separately trained NNs are applied to the respective mesh in order to approximate the last steps of the optimisation on the corresponding mesh. The result is then projected to the next finer mesh and the procedure is repeated a fixed number of times. A SIMP density field topology optimisation is directly performed in [11]. The problem can be represented in terms of the NN activation function. Different beam problems comparable to our experiments are depicted. Fully connected DNNs are used in [12] to represent implicit level set function describing the topology. For optimisation, a system of parameterised ODEs is used. A two-stage NN architecture which by construction reduces the problem of structural disconnections is developed in [13]. Deep generative models for topology optimisation problems with varying design constraints and data scenarios are explored in [14]. In [15], direct predictions without any iteration scheme and also the nonlinear elastic setting are considered. Examples are only shown for a coarse mesh discretisation of the design domain. In [16], an NN-assisted design method for topology optimisation is devised, which does not require any optimised data. A predictor NN provides the designs on the basis of boundary conditions and degree of filling as input data for which no optimisation training data are required.

The main goal of this paper is to devise new NN architectures that lower the computational burden of structural topology optimisation based on a continuous phase-field description of the density in the design domain. In particular, the approach should be able to cope with adaptive mesh refinements during the optimisation process, which has shown to significantly improve the performance of the optimisation. Moreover, as a consequence of an efficient computation in a deterministic setting, a goal is to transfer the developed techniques to the stochastic setting for the risk-averse topology optimisation task. The general strategy is to combine conventional topology optimisation methods and NNs in order to reduce the number of required iteration steps within the optimisation procedure, increasing the overall performance.

The main achievements of this paper are two new NN architectures that are demonstrated to yield state-of-the-art numerical results with a much lower number of iterations than with a classical optimisation. Moreover, in contrast to several other works that are solely founded on the image level of topology, our architectures make use of a very versatile functional phase-field description of the material distribution, which we have not observed in the literature with regard to NNs. This also holds true for the stochastic risk-averse framework, which to our knowledge has not been considered with NN predictions yet. Another novelty is the mixture of a fine reference mesh and adaptive iteration meshes during optimisation.

Inspired by the work of [5,10], as a first new NN architecture we develop a new CNN approach and show that it can replicate the reference results of [1,2]. In contrast to [10], we only have to train one NN for the entire optimisation despite mesh refinement being carried out in the iterative procedure. We subsequently extend this approach to the stochastic setting with risk-averse optimisation from [2]. Based on the CNN, we further extend the optimisation with a Long Short-Term Memory (LSTM) architecture as a second novel method. It encodes the change of the topology over several iteration steps, thus resulting in a more accurate prediction of the next gradient step.

In the numerical experiments, it can be observed that the two presented architectures perform equally well as our reference implementation. However, fewer iteration steps are required (i.e., larger steps can be used) since the gradient step predictions seem to be better than when computed with a classical optimisation algorithm.

The structure of the paper is as follows. In Section 2, we introduce the underlying setting from [1,2] and discuss the algorithms used for phase-field-based topology optimisation. In this context, we introduce the linear elasticity model and derive a state equation, the adjoint equation and a gradient equation, whose joint solution constitutes the optimisation problem under consideration. In Section 3, we present two different architectures of the NNs approximating multiple steps of the gradient equation. We start with a CNN that is well suited for processing the discretised solutions of the equations from Section 2. This is then extended to a long short-term memory NN, which is able to process a sequence of these solutions simultaneously and thus achieves a higher prediction quality. Since the data of the finite element simulation do not directly match the required structure of NNs, we provide a discussion of the data preparation for both architectures. Section 4 illustrates the practical performance of the developed NN architectures with a standard benchmark (a 2D bridge problem). The work ends with a summary and discussion of the results and some ideas for further research in Section 5. The appendix provides some background on the used problem, in particular the standard benchmark problem in Appendix A and the finite element discretisation in Appendix B. The implementation and codes for the generation of graphics and data to reproduce this work are publicly available (https://github.com/MarvinHaa/DeepNNforTopoOptisation accessed at 1 June 2022).

2. Topology Optimisation under Uncertainties

We are concerned with the task of topology optimisation with respect to a state equation of linear elasticity. This problem becomes more involved when stochastic data are assumed. In our case, this concerns material properties and the forces acting on the designed structure. These translate into the engineering world as material imperfections or fluctuations and natural forces such as wind or ocean waves. Such random phenomena are modelled in terms of random fields that often are assumed to be Gaussian with certain mean and covariance.

It is instructive to first present the deterministic topology optimisation task, which we discuss in the following Section 2.1. Subsequently, in Section 2.2 we extend the model to exhibit random data, allowing an extension of the cost functional to also include the fluctuations of the data in terms of a risk measure. In our case, this is the so-called conditional value at risk (CVaR).

For the sake of a self-contained presentation, we provide all equations that lead to the actual optimisation problem, which is given in terms of a gradient that evolves a phase field. Thus, the entire problem formulation can be understood and the required extensions to obtain the risk-averse formulation become clear. However, in case the reader is only interested in the proposed NN architectures, it might be sufficient to simply gloss over the most important parts of the problem definition, for which we provide a guideline as follows: The linear state equation is given in Equation (1), leading to the weak form in Equation (2) that is used for the computation of finite element solutions. These are required in the deterministic minimisation problem given in Equation (4), which is solved iteratively by computing the gradient step defined by Equation (6). A similar problem formulation, extended by an approximation of the CVaR risk measure, can be obtained in case of the risk-averse optimisation. This is given in Equation (9) and can be solved iteratively with gradient steps defined by Equation (11).

The presentation of this section is based on [1,2] where the optimisation problem computes the distribution of material in a given design domain described by a continuous phase field depending on the realisation of the random parameters. The optimum of this problem maximises stiffness while minimising the volume of material for the given data.

2.1. Deterministic Model Formulation

The goal is to determine an optimal distribution of a material (with density or probability) $m\in [0,1]$ in a compact design domain $D\subset {\mathbb{R}}^d,d\in \mathbb{N}$. We further assume that D satisfies sufficient regularity assumptions such that the PDE state equation exhibits a unique solution. The desired optimality of the task means that the resulting topology is as resilient (or stiff) as possible with respect to the deformation caused by the expected forces acting on it, which are described by a differentiable vector field $u:D\to {\mathbb{R}}^d$.

Definition 1.

The distribution of a material $m\in [0,1]$ in $D\subset {\mathbb{R}}^d$ is represented by a phase field $\phi :D\to \mathbb{R}$ with $0\le \phi \left(x\right)\le 1$ for all $x\in D$, where $\phi \left(x\right)=1$ if there is material at position x and $\phi \left(x\right)=0$ if there is no material at position x. At the phase transitions we allow $0<\phi \left(x\right)<1$ to ensure sufficient smoothness for phase shift. We call the evaluation of φ topology.

Note that the actual topology is reconstructed in a post-processing step by choosing some threshold in $(0,1)$ to fix the interface between material and void phase of the phase field.

2.1.1. Linear Elasticity Model

The state equation corresponding to the above problem is described by the standard linear elasticity model [17]. To define the material tensor, we first define the strain displacement (or strain tensor) $\mathcal{E}:{\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$ by

$$\mathcal{E}\left(x\right):=\frac{1}{2}(\nabla u\left(x\right)+\nabla u^T\left(x\right)),$$

which specifies the displacement of the medium in the vicinity of position x. Moreover, a so-called blend function $\omega :\mathbb{R}\to \mathbb{R}$ is given by

$$\omega \left(x\right):=max\{x^3,0\},$$

ensuring a smooth transition between the phases. According to Hooke’s law and by using the Lamé coefficients ${\mu }_{mat}>0$ and ${\lambda }_{mat}>0$, the isotropic material tensor ${\sigma }_{mat}:{\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$ for the solid phase is given by

$${\sigma }_{mat}\left(x\right):=2{\mu }_{mat}\mathcal{E}\left(x\right)+{\lambda }_{mat}\mathrm{Tr}\left(\mathcal{E}\left(x\right)\right)I.$$

This material tensor describes the acting forces between adjacent positions in the connected material, where ${\lambda }_{\mathrm{mat}}$ and ${\mu }_{\mathrm{mat}}$ are two material parameters characterising the strain–stress relationship. For the void phase, to ensure solvability of the state equation in entire domain D, we define the tensor as a fraction of the material phases. More precisely, we set ${\sigma }_{\mathrm{void}}\left(x\right):={\epsilon }^2{\sigma }_{\mathrm{mat}}\left(x\right)$ with some small $\epsilon >0$. Hence, the material tensor (or stress tensor) $\sigma :\mathbb{R}\times {\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$ is given by

$$\sigma (\phi \left(x\right),u\left(x\right)):={\sigma }_{\mathrm{mat}}\left(u\left(x\right)\right)\omega \left(\phi \left(x\right)\right)+{\sigma }_{\mathrm{void}}\left(u\left(x\right)\right)\omega (1-\phi \left(x\right)).$$

Using the material tensor $\sigma$, a force with load $g\in {\mathbb{R}}^d$ (a pressure field) and the phase field $\phi$, the displacement vector field u is described by the state equation of the standard linear elasticity model given by

$$\begin{array}{ccc}\hfill -\mathrm{div}\left[\sigma \right(\phi \left(x\right),u\left(x\right)\left)\right]& =0\hfill & \hfill \mathrm{for}\mathrm{all}x\in D,\\ \hfill u\left(x\right)& =0\hfill & \hfill \mathrm{for}\mathrm{all}x\in {\Gamma }_D,\\ \hfill \sigma \left(\phi \right(x),u(x\left)\right)& =g\hfill & \hfill \mathrm{for}\mathrm{all}x\in {\Gamma }_g,\\ \hfill u\left(x\right)·n\left(x\right)& =0\hfill & \hfill \mathrm{for}\mathrm{all}x\in {\Gamma }_s,\\ \hfill \sigma \left(\phi \right(x),u(x\left)\right)·n\left(x\right)& =0\hfill & \hfill \mathrm{for}\mathrm{all}x\in {\Gamma }_0=\partial D\setminus ({\Gamma }_D\cup {\Gamma }_g\cup {\Gamma }_s).\end{array}$$

(1)

This implies that on boundary subspace ${\Gamma }_D\subset D$ the material is fixed while on ${\Gamma }_g\subset D$ the force g acts on the material. On the boundary ${\Gamma }_s\subset D$ the material is barred from movement in normal direction n. In the following, equality is to be generally understood in a pointwise manner.

2.1.2. State Equation

The weak formulation of the state Equation (1) can be formulated as: find $u\in H_{{\Gamma }_g}^1\left(D\right)$ such that

$${\int }_D\sigma (\phi ,u)\mathcal{E}\left(v_u\right)\mathrm{d}\mu ={\int }_{{\Gamma }_g}g·v_u\mathrm{d}\mu \hspace{1em}\forall v_u\in H_0^1\left(D\right),$$

(2)

where $H^1\left(D\right)$ is the usual Sobolev space and $\mathrm{d}\mu$ the Lebesgue measure and

$$H_{{\Gamma }_g}^1\left(D\right):=\{u\in H^1\left(D\right)|\sigma (\phi ,u)=g\mathrm{on}{\Gamma }_g\}$$

and

$$H_0^1\left(D\right):=\{v_u\in H^1\left(D\right)|v_u=0\mathrm{on}{\Gamma }_0\}.$$

These definitions are in particular used for the finite element discretisation described in Appendix B.

2.1.3. Adjoint Equation

To define the optimisation problem, we introduce the Ginsburg–Landau functional $E^{\epsilon }:\mathbb{R}\to \mathbb{R}$, which serves as a penalty term for undesired variations and is defined by

$$E^{\epsilon }\left(\phi \right)={\int }_D\frac{\epsilon }{2}{|\nabla \phi |}^2+\frac{1}{2\epsilon }{\psi }_0\left(\phi \right)\mathrm{d}\mu ,$$

where $|·|$ is the Euclidean norm. This ensures that the solution to the optimisation problem can be interpreted as an actual smooth shape. The double well functional ${\psi }_0:{\mathbb{R}}^d\to \mathbb{R}$ with ${\psi }_0\left(x\right)={(\phi \left(x\right)-\phi {\left(x\right)}^2)}^2$ penalises values of $\phi$ that differ from 0 or 1 and the leading term limits the changes of $\phi$. This results in the cost functional $J^{\epsilon }:{\mathbb{R}}^d\to \mathbb{R}$ to be minimised,

$$J^{\epsilon }(\phi ,u)={\int }_{{\Gamma }_g}g·u\mathrm{d}\mu +\gamma E^{\epsilon }\left(\phi \right),\hspace{1em}\gamma >0.$$

(3)

The adaptivity parameter $\gamma$ controls the weight of the interface penalty and hence has a direct influence on the minimum respective to the characteristics of the resulting shape of $\phi$. In fact, $\gamma$ is chosen adaptively to avoid non-physical or highly porous topologies, see [2]. Additionally, we require the volume constraint ${\int }_D\phi \mathrm{d}\mu =m\left|D\right|$ with $m\in [0,1]$ to limit the amount of overall material.

The (displacement) state u from Equation (3) is obtained by solving the state Equation (2), which is used in the optimisation problem

$$\mathrm{minimize}{J}^{\epsilon }(\phi ,u)\mathrm{over}\phi \in H^1\left(D\right)$$

(4)

$$\mathrm{s}.\mathrm{t}.\mathrm{Equation}(2)\mathrm{holds},0\le \phi \left(x\right)\le 1\mathrm{for}\mathrm{all}x\in D\mathrm{and}{\int }_D\phi \left(x\right)\mathrm{d}\mu =m\left|D\right|.$$

The Allen–Cahn gradient flow approach is used to determine the solution $\phi$ for which the adjoint problem of Equation (4) is used to avoid the otherwise more costly calculation. It is shown in [2] that for $J^{\epsilon }$ the corresponding adjoint problem can be formulated as: find $p\in H^1\left(D\right)$ such that

$${\int }_D\sigma (\phi ,p))\mathcal{E}\left(v_p\right)\mathrm{d}\mu ={\int }_{{\Gamma }_g}g·v_p\mathrm{d}\mu \hspace{1em}\forall v_p\in H_0^1\left(D\right),$$

(5)

which is identical to the state equation. Hence, the respective adjoint solution p is equal to the solution u of Equation (2) and no additional system has to be solved.

2.1.4. Gradient Equation

With the solutions u respective to p one gradient step with adaptive step size $\tau$ can be characterised by the unique solution $(\phi ,\lambda )\in H^1\left(D\right)\times \mathbb{R}$ such that, for all $(v_{\phi },v_{\lambda })\in H_0^1\left(D\right)\times \mathbb{R}$,

$$\begin{array}{c}\frac{\epsilon }{\tau }{\int }_D({\phi }^{*}-{\phi }_n)v_{\phi }\mathrm{d}\mu +\epsilon \gamma {\int }_D\nabla {\phi }^{*}·\nabla v_{\phi }\mathrm{d}\mu +\frac{\gamma }{\epsilon }{\int }_D\frac{\partial }{\partial \phi }{\psi }_0\left({\phi }_n\right)v_{\phi }\mathrm{d}\mu \hfill \\ \hfill -{\int }_D\frac{\partial }{\partial \phi }\sigma (p,{\phi }_n)v_{\phi }\mathcal{E}\left(u\right)\mathrm{d}\mu +{\int }_D\lambda v_{\phi }\mathrm{d}\mu +{\int }_D({\phi }^{*}-m)v_{\lambda }\mathrm{d}\mu =0.\end{array}$$

(6)

The restriction on $0\le \phi \le 1$ for all $x\in D$ is realised by $\phi \left(x\right):=min\{max\{0,{\phi }^{*}\left(x\right)\},1\}$ in every iteration step. For the calculation of the minimum of Equation (4), the state Equation (1), the adjoint Equation (5) and subsequently the gradient Equation (6) are solved iteratively until $\phi$ converges. We always assume that solutions u and $\phi$ exist, which in fact can be observed numerically. The proposed procedure is described by Algorithm A1 where the solution of the integral equations takes place on a discretisation of D. The algorithm solves the state, adjoint and gradient equations in a loop until the solutions of the gradient equations only change slightly. The discretisation mesh is subsequently refined and the iterative process is restarted on this adjusted discretisation.

2.2. Stochastic Model Formulation

In the stochastic setting, the Lamé coefficients (determining the material properties) ${\lambda }_{\mathrm{mat}}:\Omega \to {\mathbb{R}}^{+}$ and ${\mu }_{\mathrm{mat}}:\Omega \to {\mathbb{R}}^{+}$ and the load scenarios $g:\Omega \to {\mathbb{R}}^d$ are treated as random variables on some probability space $(\Omega ,P)$. The randomness of the data is inherited by the solution of the state equation as well as the adjoint equation. As a result, the gradient step can be considered as a random distribution, see again [1,2]. The goal is to minimise the functional for the expected value of $\phi$ as well as for particularly unlikely events. For the formulation of an adequate risk-averse cost functional, we introduce the conditional value at risk (CVaR). The CVaR, a common quantity in financial mathematics, is defined for a random variable X by

$${\mathrm{CVaR}}_{\beta }\left[X\right]:=\mathbb{E}\left[X{\mathbb{1}}_{\{X>{\mathrm{VaR}}_{\beta }\left[X\right]\}}\right],$$

with ${\mathrm{VaR}}_{\beta }\left[X\right]:=inf\{t\in \mathbb{R}|P(X\le t)\ge \beta \}$ and $1>\beta \ge 0$. It characterises the expectation of the $\beta$-tail quantile distribution of X, hence accounting for bad outliers that may occur with low probability. The stochastic state equation can be formulated analogously to Equation (5) in the deterministic setting.

2.2.1. Adjoint Equation

For the risk-aware version of Equation (3) with respect to the CVaR parameter $\beta$, we define the cost $J_{\beta }^{\epsilon }\left(\phi \right):{\mathbb{R}}^d\to \mathbb{R}$ by

$$J_{\beta }^{\epsilon }\left(\phi \right)={\mathrm{CVaR}}_{\beta }\left[{\int }_{{\Gamma }_g}g·u\mathrm{d}\mu \right]+\gamma E^{\epsilon }\left(\phi \right),\hspace{1em}\gamma >0.$$

(7)

In the special case $\beta =0$, the CVaR is nothing else than the mean, i.e.,

$$J_0^{\epsilon }\left(\phi \right)=\mathbb{E}\left[{\int }_{{\Gamma }_g}g·u\mathrm{d}\mu \right]+\gamma E^{\epsilon }\left(\phi \right).$$

This results in the stochastic minimisation problem analogous to Equation (4) given by

$$\mathrm{minimise}{J}_{\beta }^{\epsilon }\left(\phi \right)\mathrm{over}\phi \in H^1\left(D\right)$$

(8)

$$\mathrm{s}.\mathrm{t}.\mathrm{Equation}(2)\mathrm{holds}\mathrm{a}.\mathrm{s}.,0\le \phi \left(x\right)\le 1\mathrm{for}\mathrm{all}x\in D\mathrm{and}{\int }_D\phi \mathrm{d}\mu =m\left|D\right|.$$

Following [2], the CVaR can be approximated in terms of the plus function. The solution of Equation (8) can hence be rewritten as

$$\underset{\phi \in H^1\left(D\right)}{min}{J}_{\beta }^{\epsilon }\left(\phi \right)=\underset{\phi \in H^1\left(D\right),t\ge 0}{min}\left(t+\frac{1}{1-\beta }\mathbb{E}\left[{\left({\int }_{{\Gamma }_g}g·u\mathrm{d}\mu \right)}_{+}\right]+\gamma E^{\epsilon }\left(\phi \right)\right).$$

(9)

An obvious approach to solve this optimisation problem is a Monte Carlo simulation, i.e., for each iteration step $n\in \mathbb{N}$ with evaluation of $u_n$, state Equation (1) have to be solved for different parameter realisations. The associated adjoint problem to Equation (9) reads

$${\int }_D\sigma (\phi ,p)\mathcal{E}\left(v_p\right)\mathrm{d}\mu =\left\{\begin{array}{cc}0,\hfill & \mathrm{if}{\int }_{{\Gamma }_g}gu\mathrm{d}\mu -t\le 0\hfill \\ {\int }_{{\Gamma }_g}{(1-\beta )}^{-1}gu\mathrm{d}\mu ,\hfill & \mathrm{else}\hfill \end{array}\right.\mathrm{a}.\mathrm{s}.$$

(10)

Consequently, the solution of Equation (10) is given by

$$p=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}{\int }_{{\Gamma }_g}gu\mathrm{d}\mu -t\le 0\hfill \\ {(1-\beta )}^{-1}u,\hfill & \mathrm{else}\hfill \end{array}\right.\mathrm{a}.\mathrm{s}.$$

2.2.2. Gradient Equation

Analogous to the deterministic approach, the gradient can be defined corresponding to Equation (9) by the solution $(\phi ,\lambda ,t)\in H^1\left(D\right)\times \mathbb{R}\times \mathbb{R}$, such that for all $(v_{\phi },v_{\lambda },v_t)\in H_0^1\left(D\right)\times \mathbb{R}\times \mathbb{R}$ the following equation holds,

$$\begin{array}{c}0=\frac{\epsilon }{{\tau }_{\phi }}{\int }_D({\phi }^{*}-{\phi }_n)v_{\phi }\mathrm{d}\mu +\frac{\epsilon }{{\tau }_t}{\int }_D(t-t_n)v_t\mathrm{d}\mu +{\int }_D\lambda v_{\phi }\mathrm{d}\mu +{\int }_D({\phi }^{*}-m)v_{\lambda }\mathrm{d}\mu \hfill \\ +\epsilon \gamma {\int }_D\nabla {\phi }^{*}·\nabla v_{\phi }\mathrm{d}\mu +\frac{\gamma }{\epsilon }{\int }_D\frac{\partial }{\partial \phi }{\psi }_0\left({\phi }_n\right)v_{\phi }\mathrm{d}\mu -{\int }_D\frac{\partial }{\partial \phi }\sigma (p,{\phi }_n)v_{\phi }\mathcal{E}\left(u\right)\mathrm{d}\mu \\ \hfill +\left\{\begin{array}{cc}{\int }_D{v}_t\mathrm{d}\mu ,\hfill & \mathrm{if}{\int }_{{\Gamma }_g}gu\mathrm{d}\mu -t\le 0\hfill \\ {\int }_D(1-\frac{1}{1-\beta })v_t\mathrm{d}\mu ,\hfill & \mathrm{else}\hfill \end{array}\right.a.s.\end{array}$$

(11)

The actual solution for one gradient step of the minimisation problem in Equation (8) with respect to (9) follows from

$$\phi \approx \frac{1}{S}\sum _{i=1}^S{\phi }_i^{*},$$

(12)

where $S\in \mathbb{N}$ is the number of samples of the Monte Carlo simulation. The larger $\beta$ chosen, the larger t becomes, and thus the number of evaluations of u for which ${\int }_{{\Gamma }_g}gu\mathrm{d}\mu -t\le 0$ holds true increases. In order to ensure a valid simulation of Equation (12), N must be chosen sufficiently large so that an adequate number of evaluations of u in each gradient step fulfils condition ${\int }_{{\Gamma }_g}gu\mathrm{d}\mu -t>0$. The described procedure is depicted in Algorithm A2 where the optimisation process is basically the same as in Algorithm A1. The central difference is that $N\in \mathbb{N}$ realisations of the optimisation problem have to be computed in each iteration. In practice, these are solved in parallel for N different $\omega \in \Omega$ and the results are then averaged. The large computational efforts caused by the slow Monte Carlo convergence are alleviated by the neural-network-based machine learning approaches presented in Section 3. In particular, gradient steps for arbitrary parameter realisations can be evaluated very efficiently and significantly fewer iterations (i.e., optimisation iterations) are required.

3. Neural Network Architectures

In modern scientific and engineering computing, machine learning techniques have become indispensable in recent years. The central goal of this work is to devise neural network architectures to facilitate an efficient computation of the risk-averse stochastic topology optimisation task. In this section, we develop two such architectures. The first one described in Section 3.1 is based on the popular convolutional neural networks (CNN) that were originally designed for the treatment of image data. In Section 3.2, a classical long short-term memory architecture (LSTM) is adapted to predict the gradient step.

The usage of deep neural networks with topology optimisation tasks have been previously examined in [4,10]. However, in contrast to other approaches, our architecture aims at a single NN that can be trained to handle arbitrarily fine meshes in terms of the requirements warranted during the topology optimisation process. More precisely, we seek an NN that predicts the gradient step ${\phi }_{n+k}\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ from Equation (6) discretised on an arbitrarily fine mesh $T_m$ at an arbitrary iteration step $n\in \mathbb{N}$ with given $k\in \mathbb{N}$, for ${\mathcal{N}}^k:{\mathbb{R}}^{1+d\times \left|V\right({\mathcal{T}}_m\left)\right|}\to {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ such that

$${\mathcal{N}}^k\left([{\phi }_n,u_{n+1}]\right)={\phi }_{n+k}.$$

(13)

Thus, the total number of iterations required for the topology optimisation iteration should ideally be reduced, resulting in improved practical performance. For the sake of a convenient presentation, we consider all other coefficients of Equation (6) as constant in the following analyses. Alternatively, one would have to increase the complexity in the number of degrees of freedom which are the weights describing the NN as well as the required training data. It can be assumed that with more information in the form of coefficients provided to the NN during training the accuracy of the resulting approximation of ${\phi }_n$ increases. Within the optimisation procedure, the actual calculation of the gradient step given in Equation (6) is performed on the basis of variable coefficients (e.g., $\tau$ and $\gamma$).

Since the discretisation ${\phi }_n\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ can be rewritten rather easily in tensor form, which represents the input of a CNN, this is the first architecture we consider in the next section.

3.1. Topology Convolutional Neural Networks (TCNN)

When using a visual representation of topologies as images (as they can be generated as output of a finite element simulation), the solution of Equation (6) can be easily transferred to the data structures used in CNNs. Consequently, predicting the gradient step with a CNN can be understood as a projection of the optimisation problem into a pixel-structured image classification problem. Here we assume that the calculation of the learned gradient step is encoded in the weights that characterise the NN.

In principle, the structure of a classical CNN consists of one or more convolutional layers followed by a pooling layer. This basic processing unit can repeat itself as often as desired. If there are at least three repetitions, we speak of a deep CNN and a deep learning architecture. In the convolutional layers a convolution matrix is applied to the input. The pooling layers are to be understood as a dimensional reduction of their input. Although common for image classification tasks, pooling layers are not used in the presented architecture.

3.1.1. TCNN Architecture

We follow the presentation of the pytorch documentation [18]. The input of a layer of the CNN architecture is a tensor $\mathcal{I}\in {\mathbb{R}}^{S\times C_{\mathrm{in}}\times H\times W}$. Here, $S\in \mathbb{N}$ is the number of input samples, in our case the evaluations of $\phi$ and u as presented in Section 3.1.2. It is therefore possible to calculate the gradient step ${\phi }_n$ from Equation (6) for several different loads g simultaneously. This way, Monte Carlo estimates become very efficient. $C_{in}\in \mathbb{N}$ corresponds to the number of input channels and each channel represents one dimension of an input ($\phi$ or u). $H\in \mathbb{N}$ and $W\in \mathbb{N}$ provide information about the dimension of the discretisation of the space D. The output of one CNN layer is specified by $\mathcal{O}\in {\mathbb{R}}^{S\times C_{out}\times H\times W}$. For fixed $s\le S$ and $i\le C_{\mathrm{out}}\in \mathbb{N}$ with $C_{\mathrm{out}}$ the number of output channels is given as

$${\mathcal{O}}_{s,i}\left({\mathcal{I}}_s\right)=b_i+\sum _{k=1}^{C_{in}}{\mathcal{W}}_{i,k}\ast {\mathcal{I}}_{s,k}.$$

(14)

Here, * denotes the cross-correlation operator, $b\in {\mathbb{R}}^{C_{out}\times H\times W}$ and ${\mathcal{I}}_{s,k}$ with $s\le S,k\le C_{\mathrm{in}}$ is a cutout of $\mathcal{I}$. The weight tensor $\mathcal{W}\in {\mathbb{R}}^{C_{\mathrm{out}}\times C_{in}\times H_K\times W_K}$ determines the dimensions of the kernel (or convolution matrix) of the layers with $H_K,W_K\in \mathbb{N}$.

For simplicity, we henceforth assume $S=1$ unless otherwise specified. In particular, the entries of the weight tensor $\mathcal{W}$ are parameters that are optimised during the training of the CNN. Depending on the architecture of the CNN, an activation function $\sigma :\mathbb{R}\to \mathbb{R}$ evaluated element-wise can additionally be applied to Equation (14).

Definition 2.

Let $b\in {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}$ and $\mathcal{W}\in {\mathbb{R}}^{C_{\mathrm{in}}\times C_{\mathrm{out}}\times H_K\times W_K}$ with $L,H,W,H_K,W_K,$$C_{\mathrm{in}},C_{\mathrm{out}}\in \mathbb{N}$ be given by one parameter vector $\theta \in {\mathbb{R}}^d$ with

$$d=C_{\mathrm{out}}·H·W+C_{\mathrm{out}}·C_{\mathrm{in}}·H_K·W_K.$$

Furthermore, let σ be a continuously differentiable activation function. We call a function

$$\mathit{Conv}(·;\theta ):{\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}\to {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}$$

(15)

a convolution layer with activation function σ if it satisfies

$$\mathit{Conv}{(\mathcal{I};\theta )}_i=\sigma \left(b_i+\sum _{k=1}^{C_{in}}{\mathcal{W}}_{i,k}\ast {\mathcal{I}}_k\right)$$

(16)

with $i=1,\dots ,C_{\mathrm{out}}$.

A sequential coupling of this layer structure provides the framework for the CNN. Specifically, for $i^L=1,\dots ,C_{out}^L\in \mathbb{N}$,

$$\begin{array}{c}\mathrm{Conv}{(·;{\theta }^{L-1})}_{i^L}^L\circ \mathrm{Conv}{(·;{\theta }^{L-1})}^{L-1}\circ \dots \circ \mathrm{Conv}{(\mathcal{I};{\theta }^1)}^1=\hfill \\ {\sigma }^L(b_{c_{\mathrm{out}}}^L+\sum _{k_L=1}^{C_{\mathrm{out}}^{L-1}}{\mathcal{W}}_{c_{\mathrm{out}},k_L}^L\ast {\sigma }^{L-1}(b_{c_{\mathrm{out}}}^{L-1}+\sum _{k_{L-1}=1}^{C_{\mathrm{out}}^{L-2}}{\mathcal{W}}_{c_{\mathrm{out}},k_{L-1}}^{L-1}\ast \dots \\ \hfill \dots \ast {\sigma }^1\left(b_{c_{\mathrm{out}}}^1+\sum _{k_1=1}^{C_{\mathrm{in}}}{\mathcal{W}}_{c_{\mathrm{out}},k_1}^1\ast {\mathcal{I}}_{k_1}\right)\dots ){)}_{i^L}.\end{array}$$

(17)

Definition 3

(CNN architecture). Let ${\mathcal{W}}^1\in {\mathbb{R}}^{C_{\mathrm{in}}\times C_{\mathrm{out}}^1\times H_K\times W_K}$, ${\mathcal{W}}^l\in {\mathbb{R}}^{C_{\mathrm{out}}^{l-1}\times C_{\mathrm{out}}^l\times H_K\times W_K}$ for $1<l\le L$ and $b^l\in {\mathbb{R}}^{C_{\mathrm{out}}^l\times H\times W}$ for $l\le L$ with $L,H,W,H_K,W_K,$$C_{\mathrm{in}},C_{\mathrm{out}}^1,\dots ,C_{\mathrm{out}}^L\in \mathbb{N}$ be given by some parameter vector $\theta \in {\mathbb{R}}^d$ with

$$d=\underset{\mathrm{Dimension}\mathrm{of}{\mathcal{W}}^1\mathrm{and}{b}^1}{\underbrace{C_{\mathrm{out}}^1\ast (C_{\mathrm{in}}(H_K·W_K)+(H·W))}}+\underset{\mathrm{Dimension}\mathrm{of}{\mathcal{W}}^l\mathrm{and}{b}^l\mathrm{for}2\le l\le L}{\underbrace{\sum _{l=2}^L{C}_{\mathrm{out}}^l\ast (C_{\mathrm{out}}^{l-1}(H_K·W_K)+(H·W))}}.$$

Furthermore, let ${\sigma }^1,\dots ,{\sigma }^L$ be given continuously differentiable activation functions. We call an NN of the form of Equation (17) an L-layer topology convolutional neural network (TCNN) and characterise it as the mapping

$${\mathcal{N}}_{\mathrm{CNN}}(·;\theta ):{\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}\to {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}.$$

The approximation of ${\mathcal{N}}_{\mathrm{CNN}}$ is hence determined by its parameter vector $\theta$. For general CNNs, the dimension $H\times W$ does not have to be constant across the different layers. The same holds true for the dimensions $H_K\times W_K$ of the kernel matrices. In fact, before implementing the convolution, we embed each channel of our input in a $(H+\lfloor \frac{H_K}{2}\rfloor )\times (W+\lfloor \frac{W_K}{2}\rfloor )$ space to preserve the dimension in the output.

Example 1.

The following specific TCNN has proved to be the most suitable for integration into Algorithm A1 for the selections of hyperparameters we have investigated. The architecture is given as an $L=6$ layer TCNN with $C_{\mathit{in}}=3$ input channels, $C_{\mathit{out}}^l=15$ for $1<l<5$ hidden channel, $C_{\mathit{out}}^6=1$ output channel and kernel size $H_K=W_K=3$ as well as trained weights described by $\theta \in {\mathbb{R}}^{8806}$ which determine the mapping by

$${\mathcal{N}}_{\mathit{CNN}}(·;\theta ):{\mathbb{R}}^{3\times 201\times 101}\to {\mathbb{R}}^{1\times 201\times 101},$$

(18)

with activation function ${\sigma }^6\left(x\right):=min\{max\{x,0\},1\}$. In contrast to many standard architectures, only the activation function of the output layer is not the identity. We chose ${\mathbb{R}}^{3\times 201\times 101}$ as input space in anticipation of the setting in Example 2, reflecting our mesh choice to discretise domain $D=[-1,1]\times [0,1]$ with $201\times 101$ nodes, which for first-order finite elements then is the dimension of the discrete functions u and φ. The architecture is depicted in Figure 1.

3.1.2. Data Preparation

On the algorithmic level, our goal is to replace the computationally costly lines 5 and 6 of all $c_m\in \mathbb{N}$ loop iterations of Algorithm A1 with a TCNN. This is not directly possible (at least for a TCNN) since the input space ${\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}$ of a TCNN does not match the mesh ${\mathcal{T}}_m$ on which the finite element discretisation and thus the optimisation of ${\phi }_n$ takes place in the current optimisation step t. Hence, it is necessary to project the evaluation of $\phi$ onto the format of a CNN. For this we define a transformation between ${\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{in}}}$ and the input tensor ${\mathbb{R}}^{H\times W\times C_{in}}$ of ${\mathcal{N}}_{\mathrm{CNN}}$. As described in Appendix B, we do not assume that the mesh ${\mathcal{T}}_m$ stays fixed in the optimisation algorithm and we instead generate a sequence of different meshes ${\mathcal{T}}_m$ by some adaptive mesh refinement, which has led to significant efficiency improvements in [1]. To obtain unique transformations between the discretisation finite element space and the input space of the NN, one can interpolate the current solutions of ${\phi }_n$ and $u_{n+1}$ from ${\mathcal{T}}_m$ onto a constant reference mesh ${\mathcal{T}}_{\mathrm{const}}$ by polynomial interpolations

$$\begin{array}{cc}\hfill p& :{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{in}}}\to {\mathbb{R}}^{H·W\times C_{\mathrm{in}}},\hfill \\ \hfill q& :{\mathbb{R}}^{H·W\times C_{\mathrm{out}}}\to {\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{out}}}.\hfill \end{array}$$

Hence, during the optimisation, the current solutions are interpolated via the operator p to the reference mesh, rendering the prediction independent from the actual adaptive mesh. After the ${\mathcal{N}}_{\mathrm{CNN}}$ prediction of the gradient step on the reference mesh, it is mapped back to the actual computation mesh via q.

Consequently, we define the reference mesh ${\mathcal{T}}_{\mathrm{const}}=(V\left({\mathcal{T}}_{\mathrm{const}}\right),E\left({\mathcal{T}}_{\mathrm{const}}\right))$ with vertices V and edges E as a graph such that $\left|V\right({\mathcal{T}}_{\mathrm{const}}\left)\right|=H·W$. Each node $v_i\in V\left({\mathcal{T}}_{\mathrm{const}}\right)$ corresponds to the values of ${\phi }_n^i={\phi }_n\left(v_i\right)\in \mathbb{R}$ and $u_n^i=u_n\left(v_i\right)\in {\mathbb{R}}^d,i\le H·W=\left|V\left({\mathcal{T}}_{\mathrm{const}}\right)\right|$ at node $v_i$. The features of the nodes can hence be interpreted as rows of a feature matrix,

$$\widetilde{{\mathcal{I}}_n}=\left[\begin{array}{cc}{\phi }_1& u_1\\ ⋮& ⋮\\ {\phi }_{n_{H·W}}& u_{n_{H·W}}\end{array}\right]\in {\mathbb{R}}^{H·W\times C_{in}}.$$

(19)

The structure of ${\mathcal{T}}_{\mathrm{const}}$ is illustrated in Figure 2.

One can now define a transformation between ${\mathbb{R}}^{H·W\times C_{in}}$ and ${\mathbb{R}}^{C_{in}\times H\times W}$ by

$$\begin{array}{cc}\hfill {\Phi }_{C_{in}}& :{\mathbb{R}}^{H·W\times C_{\mathrm{in}}}\to {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W},\hfill \\ \hfill {\Phi }_{C_{\mathrm{out}}}& :{\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}\to {\mathbb{R}}^{H·W\times C_{\mathrm{out}}}.\hfill \end{array}$$

(20)

Hence, the approximation of the gradient step 6 of Algorithm A1 is basically a coupling of the mappings p, $\Phi$ and ${\mathcal{N}}_{CNN}$, namely

$$\begin{array}{ll}{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{in}}}\stackrel{p}{\to }{\mathbb{R}}^{H·W\times C_{\mathrm{in}}}\stackrel{{\Phi }_{C_{\mathrm{in}}}}{\to }& {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}\stackrel{{\mathcal{N}}_{\mathrm{CNN}}}{\to }\\ & \stackrel{{\mathcal{N}}_{\mathrm{CNN}}}{\to }{\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}\stackrel{{\Phi }_{C_{\mathrm{out}}}}{\to }{\mathbb{R}}^{H·W\times C_{\mathrm{out}}}\stackrel{q}{\to }{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{out}}}.\end{array}$$

(21)

Example 2

(Illustrating the TCNN). The NN given by the coupling of functions in Equation (21) with ${\mathcal{N}}_{CNN}$ given as in Example 1 can be described by

$${\mathcal{N}}_{\mathrm{CNN}}^k(·;\theta ):{\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|\times 3}\to {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|},$$

(22)

with

$$\left[\begin{array}{ccc}{\phi }_n& u_{n+1}^x& u_{n+1}^y\end{array}\right]\mapsto \left[\begin{array}{c}{\phi }_{n+k}\end{array}\right]$$

(23)

for ${\phi }_n\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ and $u_n=(u_n^x,u_n^y),u_n^x,u_n^y\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ defined on ${\mathcal{T}}_m$. Hence, this NN can be applied directly to the finite element discretisations ${\phi }_n$ and $u_n$ used in Algorithms A1 and A2.

With the TCNN from the example in Equation (23) we have extended Algorithm A1. More precisely, we have inserted an NN approximation ${\mathcal{N}}_{\mathrm{CNN}}^k({\phi }_n,u_{n+1};\theta )$ in each of the $c_m\in \mathbb{N}$ steps, which predicts k iteration steps by just one evaluation. For this, the sequence $c_m$ has to be defined in advance. We leave it to future research to adaptively control the sequence $c_m$ dynamically within the optimisation algorithm. This extension of Algorithm A1 is described by Algorithm 1. In an analogous way, we also extend Algorithm A2 by the TCNN given in Equation (23). In particular, we are able to evaluate all samples $S\in \mathbb{N}$ in parallel by adding additional sample dimensions to the input tensor of the TCNN given in Equation (14). This procedure is illustrated in Algorithm 2. Again, the parameter $c_m$ has to be chosen in advance.

Algorithm 1: Deterministic optimisation algorithm with TCNN approximated gradient step

Algorithm 2: Stochastic optimisation algorithm with TCNN approximated gradient step

3.2. Topology Long Short-Term Memory Neural Networks (TLSTM)

One possible approach to improve the prediction of ${\phi }_n$ using an NN is to provide the classifier not just one tuple $({\phi }_n,u_{n+1})$ as an inpu, but to have it process a larger amount of information by a sequence of these tuples of the last $T\in \mathbb{N}$ iteration steps, i.e.,

$$\left(({\phi }_{n-T},u_{n-T+1}),({\phi }_{n-T+1},u_{n-T+2}),\dots ,({\phi }_n,u_{n+1})\right).$$

Through this, the shift of the phase field or the change of the topology ${\phi }_n$ over time is also transferred as input to the NN. The sequence prediction problem considered in this case differs from the single step time prediction in the sense that the prediction target is now a sequence that contains both spatial and temporal information. Theoretically, this information can also be learned directly from the NN. However, in practice it is more effective to adapt the architecture to the information we have in advance (in our case with respect to the time dependency) to achieve better results. An NN that allows exactly this is a recurrent Neural Network (RNN). Unfortunately, standard RNNs often suffer from the vanishing gradient problem [19,20] which we try to prevent from the beginning. Therefore, we build on the special RNN concept of a Long Short-Term Memory (LSTM) in the context of our problem, which is more robust against the vanishing gradient issue and provides promising results, especially in the analysis of time series. For a background on time series analysis and a review of the different methods, we refer to the survey article [21]. In practice, time series are usually stored as one-dimensional sequences in vector format. Consequently, there is no out-of-the-box LSTM layer implementation for structures such as the input tensor we require in Equation (14). Nevertheless, we still do not wish to abandon the mechanism of convolution within the NN in order to keep the structural information of $\phi$ and u. An LSTM layer with a convolutional structure can be constructed by replacing the matrix vector multiplication within a standard LSTM layer by convolutional layers. The unique advantage of an LSTM according to [20] is its cell-gate architecture, which mitigates the vanishing gradient problem. More precisely, it consists of a “memory cell” $c_t:{\mathbb{R}}^{4\times d_{\mathrm{in}}^t}\to {\mathbb{R}}^{d_{\mathrm{out}}^t}$ that serves as an accumulator of the current state $t\le T,t\in \mathbb{N}$, in the processed sequence. The information capacity of the last status $c_{t-1}$ within $c_t$ is controlled by the activation of the so-called “forget gate” $f_t:{\mathbb{R}}^{3\times d_{\mathrm{in}}^t}\to {\mathbb{R}}^{d_{\mathrm{out}}^t}$. The information capacity of the input state $x_t\in {\mathbb{R}}^{d_{\mathrm{in}}^t}$ is controlled by the activation of the input gate $i_t$. Which information (or whether any at all) gets transferred from memory cell $c_t$ to state $h_t:{\mathbb{R}}^{2\times d_{\mathrm{in}}^t}\to {\mathbb{R}}^{d_{\mathrm{out}}^t}$ is in turn controlled by the activation of the output gate $o_t$. From a technical point of view, the gates can be understood as learning forward layers.

3.3. TLSTM Architecture

An ordinary LSTM layer to generate complex sequences with long-range structure as presented in [22] corresponds to the described logic above and can be formulated numerically for a sequence of one-dimensional input state $x_t\in {\mathbb{R}}^{d_{\mathrm{in}}^t}$ and output vector $h_t\in {\mathbb{R}}^{d_{\mathrm{out}}^t}$ as an equation system

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill i_t(x_t,h_{t-1},c_{t-1})& =\sigma ({\mathcal{W}}_{xi}{x}_t+{\mathcal{W}}_{hi}{h}_{t-1}+{\mathcal{W}}_{ci}\odot c_{t-1}+b_i),\hfill \\ \hfill f_t(x_t,h_{t-1},c_{t-1})& =\sigma ({\mathcal{W}}_{xf}{x}_t+{\mathcal{W}}_{hf}{h}_{t-1}+{\mathcal{W}}_{cf}\odot c_{t-1}+b_f),\hfill \\ \hfill c_t(f_t,c_{t-1},i_t,x_t,h_{t-1})& =f_t\odot c_{t-1}+i_t\odot tanh({\mathcal{W}}_{xc}{x}_t+{\mathcal{W}}_{hc}{h}_{t-1}+b_c),\hfill \\ \hfill o_t(x_t,h_{t-1},c_t)& =\sigma ({\mathcal{W}}_{xo}{x}_t+{\mathcal{W}}_{ho}{h}_{t-1}+{\mathcal{W}}_{co}\odot c_t+b_o),\hfill \\ \hfill h_t(o_t,c_t)& =o_t\odot tanh\left(c_t\right),\hfill \end{array}\end{array}$$

(24)

where $\sigma :\mathbb{R}\to \mathbb{R}$ with $\sigma \left(x\right)=\frac{1}{1+e^x}$ and tanh are evaluated element-wise. The operation ⊙ denotes the Hadamard product and the subscripts of the weight matrices $\mathcal{W}\in {\mathbb{R}}^{d_{\mathrm{out}}^t\times d_{\mathrm{in}}^t}$ describe the affiliation to the gates. For example, ${\mathcal{W}}_{xi}$ is the weight matrix to input $x_t$ of gate $i_t$. This illustrates how the weights of the LSTMs are transferred to the weights of convolution LSTMs in the following.

We want to reformulate Equation (24) by replacing all matrix-vector multiplications (i.e., the forward layer) with a convolution layer from Definition 2. This is inspired by [23], which has already provided the theoretic architecture of a convolutional LSTM layer with the approach on precipitation forecasting. Let $\mathrm{Conv}(·;\theta ):{\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}\to {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}$ be a convolutional layer and $\mathcal{I}\in {\mathbb{R}}^{T\times C_{\mathrm{in}}\times H\times W}$ a sequence of inputs ordered by the discrete time dimension $T\in \mathbb{N}$. A convolutional LSTM layer to an input sequence ${\mathcal{I}}_{t\le T}$ and ${\mathcal{H}}_0=0\in {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W},{\mathcal{C}}_0=0\in {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}$ (since at $t=1$ we do not yet have any information about earlier steps in the sequence) is given by a system of equations,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{i}}_t({\mathcal{I}}_t,{\mathcal{H}}_{t-1},{\mathcal{C}}_{t-1})& =\sigma \left(\mathrm{Conv}({\mathcal{I}}_t;{\theta }_{xi})+\mathrm{Conv}({\mathcal{H}}_{t-1};{\theta }_{hi})+{\mathcal{W}}_{ci}\odot {\mathcal{C}}_{t-1}\right),\hfill \\ \hfill \widehat{f_t}({\mathcal{I}}_t;{\mathcal{H}}_{t-1},{\mathcal{C}}_{t-1})& =\sigma \left(\mathrm{Conv}({\mathcal{I}}_t;{\theta }_{xf})+\mathrm{Conv}({\mathcal{H}}_{t-1};{\theta }_{hf})+{\mathcal{W}}_{cf}\odot {\mathcal{C}}_{t-1}\right),\hfill \\ \hfill {\mathcal{C}}_t(\widehat{f_t},{\mathcal{C}}_{t-1},{\widehat{i}}_t,{\mathcal{I}}_t)& =\widehat{f_t}\odot {\mathcal{C}}_{t-1}+{\widehat{i}}_t\odot tanh\left(\mathrm{Conv}({\mathcal{I}}_t;{\theta }_{xc})+\mathrm{Conv}({\mathcal{H}}_{t-1};{\theta }_{hc})\right),\hfill \\ \hfill \widehat{o_t}({\mathcal{I}}_t,{\mathcal{H}}_{t-1},{\mathcal{C}}_t)& =\sigma \left(\mathrm{Conv}({\mathcal{I}}_t;{\theta }_{xo})+\mathrm{Conv}({\mathcal{H}}_{t-1};{\theta }_{ho})+{\mathcal{W}}_{co}\odot {\mathcal{C}}_t\right),\hfill \\ \hfill {\mathcal{H}}_t(\widehat{o_t},{\mathcal{C}}_t)& =\widehat{o_t}\odot tanh\left({\mathcal{C}}_t\right),\hfill \end{array}\end{array}$$

(25)

with $t\le T,t\in \mathbb{N}$ and $\mathcal{H}\in {\mathbb{R}}^{T\times C_{\mathrm{out}}\times H\times W}$ the output of the convolutional LSTM layer as well as ${\mathcal{W}}_{ci},{\mathcal{W}}_{cf}\in {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W},{\mathcal{W}}_{co}\in {\mathbb{R}}^{C_{out}\times H\times W}$ in Equation (24). The subscript t indicates a cutout of the t-th element of sequence dimension T of the respective tensor.

Definition 4

(LSTM layer). Let $L,H,W,T,C_{\mathrm{in}},C_{\mathrm{out}}\in \mathbb{N}$ as well as ${\mathcal{W}}_{ci},{\mathcal{W}}_{cf}\in {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W},{\mathcal{W}}_{co}\in {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}$ and parameter vectors, specifying the convolutional layer as in Definition 3 for $L=1$ from Equation (25),

$${\theta }_{xi}\in {\mathbb{R}}^{d_{xi}},{\theta }_{hi}\in {\mathbb{R}}^{d_{hi}},{\theta }_{xc}\in {\mathbb{R}}^{d_{xc}},{\theta }_{hc}\in {\mathbb{R}}^{d_{hc}},$$

$${\theta }_{xf}\in {\mathbb{R}}^{d_{xf}},{\theta }_{hf}\in {\mathbb{R}}^{d_{hf}},{\theta }_{xo}\in {\mathbb{R}}^{d_{xo}},{\theta }_{ho}\in {\mathbb{R}}^{d_{ho}},$$

and described by the parameter vector $\theta \in {\mathbb{R}}^d$, with

$$d=d_{xi}+d_{hi}+d_{xc}+d_{hc}+d_{xf}+d_{hf}+d_{xo}+d_{ho}+2·C_{\mathrm{in}}·H·W+·C_{out}·H·W.$$

Furthermore, let $\sigma :\mathbb{R}\to \mathbb{R}$ with $\sigma \left(x\right)=\frac{1}{1+e^x}$ and tanh be evaluated element-wise. We call a function,

$$\mathit{LSTM}(·;\theta ):{\mathbb{R}}^{3\times T\times C_{\mathrm{in}}\times H\times W}\to {\mathbb{R}}^{2\times T\times C_{\mathrm{out}}\times H\times W},$$

an LSTM layer, if it satisfies the mapping rule given by the system of Equation (25).

For the forecasting of our gradient sequence, we use an encoder–decoder architecture (i.e., an “autoencoder”) consisting of $2L,L\in \mathbb{N}$, LSTM layers,

$${\mathrm{LSTM}}^l(·;{\theta }^l):{\mathbb{R}}^{3\times T\times C_{\mathrm{in}}^l\times H\times W}\to {\mathbb{R}}^{2\times T\times C_{\mathrm{out}}^l\times H\times W},$$

that satisfies the mapping rule given by the equation system (25) with $1\le l\le 2L$. Therefore, the encoding and decoding blocks of the autoencoder have the same number of layers $L\in \mathbb{N}$. The autoencoder for an input sequence $\mathcal{I}\in {\mathbb{R}}^{T\times C_{in}\times H\times W}$ can be described by the following system of equations of the encoder block,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{LSTM}}^1({\mathcal{I}}_t,{\mathcal{C}}_{t-1}^1,{\mathcal{H}}_{t-1}^1;{\theta }^1)& =[{\mathcal{C}}_t^1,{\mathcal{H}}_t^1],\hfill \\ \hfill {\mathrm{LSTM}}^2({\mathcal{H}}_t^1,{\mathcal{C}}_{t-1}^2,{\mathcal{H}}_{t-1}^2;{\theta }^2)& =[{\mathcal{C}}_t^2,{\mathcal{H}}_t^2],\hfill \\ \hfill & ⋮\hfill \\ \hfill {\mathrm{LSTM}}^{L-1}({\mathcal{H}}_t^{L-2},{\mathcal{C}}_{t-1}^{L-1},{\mathcal{H}}_{t-1}^{L-1};{\theta }^{L-1})& =[{\mathcal{C}}_t^{L-1},{\mathcal{H}}_t^{L-1}],\hfill \\ \hfill {\mathrm{LSTM}}^L({\mathcal{H}}_t^{L-1},{\mathcal{C}}_t^L,{\mathcal{H}}_{t-1}^L;{\theta }^L)& =[{\mathcal{C}}_t^L,{\mathcal{H}}_t^L],\hfill \end{array}\end{array}$$

(26)

with $1\le t\le T_{\mathrm{en}},t\in \mathbb{N}$. This is combined with the following decoder block by setting $\widetilde{{\mathcal{O}}_0}={\mathcal{H}}_{T_{\mathrm{en}}}^L$ and ${\mathcal{H}}_0^{L+1}={\mathcal{H}}_{T_{\mathrm{en}}}^1,\dots ,{\mathcal{H}}_0^{2L}={\mathcal{H}}_{T_{\mathrm{en}}}^L$ and ${\mathcal{C}}_0^{L+1}={\mathcal{C}}_{T_{\mathrm{en}}}^1,\dots ,{\mathcal{C}}_0^{2L}={\mathcal{C}}_{T_{\mathrm{en}}}^L$,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{LSTM}}^{L+1}({\tilde{\mathcal{O}}}_{t-1},{\mathcal{C}}_{t-1}^{L+1},{\mathcal{H}}_{t-1}^{L+1};{\theta }^{L+1})& =[{\mathcal{C}}_t^{L+1},{\mathcal{H}}_t^{L+1}],\hfill \\ \hfill {\mathrm{LSTM}}^{L+2}({\mathcal{H}}_t^{L+1},{\mathcal{C}}_{t-1}^{L+2},{\mathcal{H}}_{t-1}^{L+2};{\theta }^{L+2})& =[{\mathcal{C}}_t^{L+2},{\mathcal{H}}_t^{L+2}],\hfill \\ \hfill & ⋮\hfill \\ \hfill {\mathrm{LSTM}}^{2L-1}({\mathcal{H}}_t^{2L-2},{\mathcal{C}}_{t-1}^{2L-1},{\mathcal{H}}_{t-1}^{2L-1};{\theta }^{2L-1})& =[{\mathcal{C}}_t^{2L-1},{\mathcal{H}}_t^{2L-1}],\hfill \\ \hfill {\mathrm{LSTM}}^{2L}({\mathcal{H}}_t^{2L-1},{\mathcal{C}}_t^{2L},{\mathcal{H}}_{t-1}^{2L};{\theta }^{2L})& =[{\mathcal{C}}_t^{2L},{\tilde{\mathcal{O}}}_t],\hfill \end{array}\end{array}$$

(27)

with $1\le t\le T_{\mathrm{dec}},t\in \mathbb{N}$.

It should be mentioned that the input and output sequences do not have to be of the same length (in fact, in general $T_{\mathrm{en}}\ne T_{\mathrm{dec}}$). Furthermore, $C_{\mathrm{out}}^{l-1}=C_{\mathrm{in}}^l$ holds for individual LSTM layers $2\le l\le 2L$ defined in Equations (26) and (27). Especially, since the output sequence $\tilde{\mathcal{O}}$ is also input of ${\mathrm{LSTM}}^{L+1}$, it holds that $C_{\mathrm{out}}^{2L}=C_{\mathrm{in}}^{L+1}$. In order to be able to select the dimensions of the output tensor $\tilde{\mathcal{O}}$ completely independently of the hidden channels, we additionally apply a convolutional layer $\mathrm{Conv}(·;{\theta }^{\mathrm{final}}):{\mathbb{R}}^{C_{\mathrm{out}}^{2L}\times H\times W}\to {\mathbb{R}}^{C_{\mathrm{final}}\times H\times W}$ with activation function ${\sigma }_{\mathrm{final}}:\mathbb{R}\to \mathbb{R}$ to concatenate the hidden channel to an arbitrary number of output channels $C_{\mathrm{final}}\in \mathbb{N}$ given by

$$\mathrm{Conv}({\tilde{\mathcal{O}}}_t;{\theta }^{\mathrm{final}})={\mathcal{O}}_t,\hspace{1em}1\le t\le T_{\mathrm{dec}}.$$

(28)

Definition 5

(TLSTM architecture). Let $L,H,W,T_{\mathrm{en}},T_{\mathrm{dec}},C_{\mathrm{final}},C_{\mathrm{in}}^1\in \mathbb{N}$ as well as $C_{\mathrm{out}}^l\in \mathbb{N}$, with $1\le l\le 2L$ be given. Hence, the respective LSTM layers from the encoder defined in Equation (26) and decoder in (27) block as well as the output layer in Equation (28) can be described by the parameter vectors of the CNN and LSTM layers (see Definition 3 for $L=1$ and Definition 4) ${\theta }^{\mathit{final}}\in {\mathbb{R}}^{d^f},{\theta }^l\in {\mathbb{R}}^{d^l}$ with $d^f\in \mathbb{N},d^l\in \mathbb{N}$ for $1\le l\le 2L$ and an activation function ${\sigma }_{\mathit{final}}:\mathbb{R}\to \mathbb{R}$ of the output layer. These parameter vectors as well as the weight tensors ${\mathcal{W}}_{ci}^l,{\mathcal{W}}_{cf}^l\in {\mathbb{R}}^{C_{\mathrm{in}}\times H\times W}$ and ${\mathcal{W}}_{co}^l\in {\mathbb{R}}^{C_{\mathrm{out}}\times H\times W}$ for $l\le 2L$ can in turn be described collectively by the parameter vector $\theta \in {\mathbb{R}}^d$, where

$$d=d_f+\sum _{l=1}^{2L}{d}^l.$$

We call an NN as described in Equations (26)–(28) Convolutional Topology Long Short-Term Memory (TLSTM) and characterise it by

$${\mathcal{N}}_{\mathrm{LSTM}}(·;\theta ):{\mathbb{R}}^{T_{\mathrm{en}}\times C_{\mathrm{in}}\times H\times W}\to {\mathbb{R}}^{T_{\mathrm{dec}}\times C_{\mathrm{final}}\times H\times W}.$$

The difference between a TLSTM and an LSTM is therefore the structure of the input tensor ${\mathbb{R}}^{T_{\mathrm{en}}\times C_{in}\times H\times W}$ of a TLSTM instead of ${\mathbb{R}}^{T_{\mathrm{en}}\times C_{\mathrm{in}}\times H}$ and the internal calculation carried out with convolutional layers instead of standard multiplications.

Example 3.

For the experiments in Section 4.2, the underlying $L=4$ layer TLSTM with $C_{\mathit{in}}=3$ input channels, $C_{\mathit{out}}^l=9,l\le 4$ hidden channels, $C_{\mathit{final}}=1$ output channel, kernel size of all included convolutional layers $H_K=W_K=3$ and sequence lengths $T_{\mathit{en}}=5$, $T_{\mathit{dec}}=10$ with trained weights described by $\theta \in {\mathbb{R}}^{509266}$ are given as

$${\mathcal{N}}_{\mathrm{LSTM}}·;\theta ):{\mathbb{R}}^{5\times 3\times 201\times 101}\to {\mathbb{R}}^{10\times 1\times 201\times 101},$$

(29)

with activation function ${\sigma }^{\mathit{final}}\left(x\right)=min\{max\{x,0\}1\}$. The autoencoder structure for an 8-layer LSTM described in Equations (26)–(28) for (29) is visualised in Figure 3.

3.4. Data Preparation

As in the case of the integration of the ${\mathcal{N}}_{\mathrm{CNN}}$ proposed in Section 3.1, we want to replace lines 5 and 6 in Algorithm A1 with the approximation of the ${\mathcal{N}}_{\mathrm{LSTM}}$ from Definition 5. In case of a TLSTM, the evaluations of $\phi$ and u have to be transformed from the graph structure of the finite element simulation into an appropriate tensor format. This in principle is analogous to the composition ${\Phi }_{C_{\mathrm{in}}}\circ p$ in Equation (21). The only difference is that now this is performed on a sequence of evaluations ${\phi }_{m_n}$ and $u_{m_n}$ of length $T_{\mathrm{en}}\in \mathbb{N}$. As the subscripts suggest, such a sequence does not necessarily have to be evaluated on a fixed mesh ${\mathcal{T}}_m$, it may extend over a sequence of meshes ${\mathcal{T}}_m$. However, since we use polynomial interpolations

$$\begin{array}{cc}\hfill p_T& :{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{in}}\times T}\to {\mathbb{R}}^{H·W\times C_{\mathrm{in}}\times T},\hfill \\ \hfill q_T& :{\mathbb{R}}^{H·W\times C_{\mathrm{out}}\times T}\to {\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{out}}\times T}\hfill \end{array}$$

to transfer the sequence ${\phi }_{n-T},\dots ,{\phi }_n$ and $u_{n-T+1},\dots ,u_{n+1}$ onto some reference mesh ${\mathcal{T}}_{\mathrm{const}}=(V\left({\mathcal{T}}_{\mathrm{const}}\right),E\left({\mathcal{T}}_{\mathrm{const}}\right))$.

We intend to process $T_{\mathrm{en}}$ feature matrices of the form of Equation (19). Hence, we define transformations

$${\Phi }_{C_{\mathrm{in}},T}:{\mathbb{R}}^{H·W\times C_{\mathrm{in}}\times T}\to {\mathbb{R}}^{T\times C_{\mathrm{in}}\times H\times W},$$

(30)

$${\Phi }_{C_{\mathrm{out}},T}:{\mathbb{R}}^{T\times C_{\mathrm{out}}\times H\times W}\to {\mathbb{R}}^{H·W\times C_{\mathrm{out}}\times T}.$$

Thus, the approximation of a gradient step in Algorithm A1 by a TLSTM can be understood as a concatenation of the form

$$\begin{array}{ll}{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{in}}\times T_{\mathrm{en}}}& \stackrel{p_T}{\to }{\mathbb{R}}^{H·W\times C_{\mathrm{in}}\times T_{\mathrm{en}}}\stackrel{{\Phi }_{C_{\mathrm{in}},T_{\mathrm{en}}}}{\to }{\mathbb{R}}^{T\times C_{\mathrm{in}}\times H\times W}\stackrel{{\mathcal{N}}_{\mathrm{CNN}}}{\to }\\ & \stackrel{{\mathcal{N}}_{\mathrm{CNN}}}{\to }{\mathbb{R}}^{T\times C_{\mathrm{out}}\times H\times W}\stackrel{{\Phi }_{C_{\mathrm{out}},T_{\mathrm{dec}}}}{\to }{\mathbb{R}}^{H·W\times C_{\mathrm{out}}\times T_{\mathrm{dec}}}\stackrel{q_T}{\to }{\mathbb{R}}^{|V\left({\mathcal{T}}_m\right)|\times C_{\mathrm{out}}\times T_{\mathrm{dec}}}\end{array}$$

(31)

Example 4

(The TLSTM). The NN given by the coupling of functions from Equation (31), where ${\mathcal{N}}_{\mathit{LSTM}}$ as given in Example 3, can be described by

$${\mathcal{N}}_{\mathit{LSTM}}(·;\theta ):{\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|\times 3\times 5}\to {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|\times 10}$$

(32)

and

$$\left[\begin{array}{ccc}{\phi }_{n-5}& u_{n-5+1}^x& u_{n-5+1}^y\\ & ⋮& \\ {\phi }_n& u_{n+1}^x& u_{n+1}^y\end{array}\right]\mapsto \left[\begin{array}{c}{\phi }_{n+1}\\ ⋮\\ {\phi }_{n+10}\end{array}\right]$$

for ${\phi }_n\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ and $u_n=(u_n^x,u_n^y),u_n^x,u_n^y\in {\mathbb{R}}^{\left|V\right({\mathcal{T}}_m\left)\right|}$ defined on mesh ${\mathcal{T}}_m$.

The TLSTM from Example 4 can directly be integrated into Algorithm A1 as with the previous integration with Algorithm 1. In fact, since in Algorithm A1 only the most recent gradient step is relevant, in practice we restrict the inverse mapping on the last element of the predicted sequences to save calculation time. The only difference is that the sequences ${\phi }_{n-T},\dots ,{\phi }_n$ and $u_{n-T+1},\dots ,u_{n+1}$ have to be stored in a list. We chose $c_1=125$ and $c_m=50$ for all $2\le m\in \mathbb{N}$. This procedure is described in Algorithm 3. As in the case of the TCNN, we are able to include the extra sample dimension S to approximate gradient steps from multiple problems at once.

Algorithm 3: Deterministic optimisation algorithm with TLSTM approximated gradient step

4. Results

This section is devoted to numerical results of the two previously described neural network architectures. The implementations were conducted with the open source packages PyTorch [18] for the NN part and FEniCS [24] for the FE simulations (see introduction for the link to the code repository). We first illustrate the performance of the TCNN in Section 4.1 with a deterministic bridge example compared to a classical optimisation. The important observation is that with the TCNN the optimisation can be carried out with far fewer optimisation steps while still leading to the reference topologies from [1]. Similar results can be observed for the risk-averse stochastic optimisation. In Section 4.2, numerical experiments of the TLSTM architecture are presented. The performance is revealed to be comparable to the TCNN architecture and the optimisation appears to be more robust with respect to the data realisations.

4.1. TCNN Examples

Before we can use the TCNN architecure for the optimisation in Algorithm 1, we have to train it on data that describe the system response of Equation (6). Note that it would not be useful to allow the ${\mathcal{N}}_{\mathrm{CNN}}$ to learn the gradient steps of a fixed setting since different settings of the bridge problem from Appendix A.1 should efficiently be tackled. In considering the stochastic setting of the problem as defined in Section 2.2, the TCNN is trained to learn the gradient steps $\phi$ for a random $g:\Omega \to {\mathbb{R}}^d$.

4.1.1. Sampling the Data

To train the architecture, appropriate training data have to be generated. In order to achieve this, we chose the same setting for g as in Appendix A.2. Using the optimiser in Algorithm A1, we can generate $S\in \mathbb{N}$ different sample paths of gradient steps ${\phi }_n\left(g\right)$ and solutions of the state equation $u_n\left(g\right)$ by generating S samples of g. In this procedure, we store every $k\in \mathbb{N}$ iteration step of ${\phi }_n$ and $u_n$ in order to approximate k gradient steps at once. More precisely, we store $\lfloor \frac{N_{\max }}{k}\rfloor -1$ tuples

$$\left(\left[\begin{array}{ccc}{\phi }_n& u_{n+1}^x& u_{n+1}^y\end{array}\right],\left[\begin{array}{c}{\phi }_{n+k}\end{array}\right]\right),$$

(33)

with $0\le n\le N_{\max }-k$, where $N_{\max }\in \mathbb{N}$ is fixed in advance, representing the maximum number of iterations of an optimisation. For the training of the models in the following experiments, we have chosen $N_{\max }=500$ since the topologies have mostly converged after this number of iterations. The overall number of $S(\lfloor \frac{N_{\max }}{k}\rfloor -1)$ tuples are merged into an unsorted data set ${\mathcal{D}}_{\mathrm{CNN}}$.

4.1.2. TCNN Predictions

The following experiment validates that the performance of Algorithm A1 can be replicated or (desirably) improved by including a CNN as described in Algorithm 1. As a first test, we illustrate that the proposed new architecture is indeed capable of predicting the gradients of the optimisation procedure.

Figure 4 shows the evaluations of the model Equation (22) after determining $\theta$ within the training of the NN on the data set ${\mathcal{D}}_{CNN}$. Here, the prediction ${\mathcal{N}}_{\mathrm{CNN}}^k({\phi }_n,u_{n+1};\theta )$ and the actual gradient step ${\phi }_{n+k}$ generated by Algorithm A1 are compared for different loads sampled from a truncated normally distributed g. Since the predictions of ${\mathcal{N}}_{\mathrm{CNN}}^5$ are hardly distinguishable from the reference fields, we have also trained Equation (22) to predict larger time steps (for 25 and 100 iteration steps at once). However, these NNs have proved to be less reliable in practice as prediction quality decreases. An illustrative selection of some predictions is provided in Figure A5 and Figure A6 in Appendix C.

4.1.3. Deterministic Bridge Optimisation

In these experiments we compare the performance of Algorithm 1 with that of Algorithm A1 in the setting of Appendix A.1. For the NN in Algorithm 1, we use Equation (22) from Example 2. For $c_m\in \mathbb{N}$ in Algorithm 1 we chose $c_m=55$ for all $m\in \mathbb{N}$. In order to train Equation (22), data of the form of Equation (33) from Section 4.1.1 are used. We expect that the best (reference) results in the setting from experiment A.1 can be obtained, since the distribution of the training data is a truncated normal distribution around this expected value and we thereby have an accumulation of training points around the load $g={(0,-5000)}^T$. As a convergence criterion, we used the convergence criterion of the mesh refinement of Equation (A1) on a maximally fine mesh with $\left|V\right({\mathcal{T}}_m\left)\right|\le 15,000$. A sub-sequence generated by Algorithm 1 is shown in Figure A7a in comparison to that of Appendix A.1 in Figure A7b in Appendix C. There, it can be observed that the classical and CNN-assisted optimisation results essentially appear identical with the faster CNN convergence.

The metrics for evaluating our algorithms have also improved through the application of the CNN as can be observed in Figure 5. For better comparability, we ran both algorithms 10 times and averaged all metrics. To be more precise, this is only the average of the calculation time, as the remaining metrics are deterministic and therefore always the same. It is easy to see how the metrics diverge with the first application of the CNN at iteration step 55, especially by the computation time required per iteration. The most significant indicator is the evaluation of $J^{\epsilon }\left({\phi }_n\right)$ per iteration of Equation (3) at the top-right of Figure 5. The graph for Algorithm 1 reaches a constant lower level than that of Algorithm A1 after about 250 iterations and thus fulfils the convergence criterion earlier. Accordingly, the step size criterion for ${\tau }_n$ applies earlier by using the CNN, which further accelerates convergence. An interesting insight is provided by the calculation time, which shows that the actual time required per iteration step is more or less the same, except for the iteration steps in which Equation (22) is applied. This is indicated by the upward outliers in the computation time series. This additional computation cost can be explained by the application of the mesh projection of Equation (21), which represents an aspect requiring further improvements. Nevertheless, in total we achieve a shorter total run-time due to the faster convergence of Algorithm 1.

A detailed list of the run-times and target value metrics for the functional $J^{\epsilon }\left({\phi }_n\right)$ is provided in

Table 1. Comparison of metrics between Algorithms A1 and 1.

Method	Applied Load g	Converges after ${\mathit{n}}_{\mathbf{final}}$	Evaluation of ${\mathit{J}}^{\mathit{\epsilon }}\left({\mathit{\phi }}_{{\mathit{n}}_{\mathbf{final}}}\right)$	Compliance ${\int }_{{\mathbf{\Gamma }}_{\mathit{g}}}\mathit{g}·\mathit{u}\left(\mathit{\phi }\right)\mathbf{d}\mathit{s}$
Algorithm A1	${(0,-5000)}^T$	538	1475.39	1173.22
Algorithm 1	${(0,-5000)}^T$	441	1206.49	1148.66

. The evaluation of $J^{\epsilon }\left({\phi }_{n_{\mathrm{final}}}\right)$ denotes the value of the functional for the topology ${\phi }_{n_{\mathrm{final}}}$ converged after $n_{\mathrm{final}}\in \mathbb{N}$ iteration steps. The compliance is the value that is actually minimised in terms of the functional $J^{\epsilon }\left({\phi }_{n_{\mathrm{final}}}\right)$. Algorithm 1 requires less computing time than the reference procedure after extending it with the CNN architecture from Equation (22).

4.1.4. Stochastic Bridge Optimisation

Algorithm A2 can easily be extended by the TCNN of Equation (22) in order to improve the efficiency for topology optimisation under uncertainties. The corresponding procedure is shown in Algorithm 2 where we chose $c_m=55$ for all $m\in \mathbb{N}$. Note that the predictions of the different realisations of ${\phi }_n\left({\omega }_i\right),{\omega }_i\in \Omega$ for $i=1,\dots ,S\in \mathbb{N}$ (where an evaluation ${\phi }_n\left({\omega }_i\right)$ are to be interpreted as transformation of an evaluation from $g\left({\omega }_i\right)$) are actually not executed within a loop but in parallel (lines 8–12 of Algorithm 2). This is possible because NNs are generally able to process batches of data in parallel. We have also implemented parallelisation for Algorithm A2, which is limited by the number of processor cores of the actual compute cluster. We want to compare the performance of Algorithm A2 and Algorithm 2 using the same setting as in Appendix A.2. To ensure comparable results despite stochastic parameters, we set the random seed to 42 before running both algorithms. The resulting sub-sequences of ${\phi }_n$ are compared side by side in Figure 6. Although the topology converges after fewer iterations with Algorithm A2, one can see that the topology resulting from Algorithm 2 has a more stable shape since the topology does not lose material to the unnecessary extra spoke. This is confirmed by the metrics in

Table 2. Comparison of metrics between Algorithms A2 and 2.

Method	Converges after ${\mathit{n}}_{\mathbf{final}}$	Evaluation of ${\mathit{J}}_0^{\mathit{\epsilon }}\left({\mathit{\phi }}_{{\mathit{n}}_{\mathbf{final}}}\right)$	Compliance $\mathbb{E}\left[{\int }_{{\mathbf{\Gamma }}_{\mathit{g}}}\mathit{g}·\mathit{u}\left({\mathit{\phi }}_{{\mathit{n}}_{\mathbf{final}}}\right)\mathbf{d}\mathit{s}\right]$	Samples
Algorithm A2	490	765.85	729.20	224
Algorithm 2	517	744.11	708.61	224
Algorithm 2	500	760.26	723.24	224

where one can see that Algorithm 2 achieved a lower compliance after fewer iteration steps. Additionally, the optimisation of Algorithm 2 is stopped after 500 iterations to show that it achieves a better result in less time as shown in Figure 7. A notable observation is that the times of applying Equation (22) in Algorithm 2 can be identified by the spikes in the computation time of the iterations. It can be seen that despite the additional time required by the transformation in Equation (21), the calculation of a stochastic gradient step using Equation (22) is generally faster. This is due to the dynamic parallelisation that PyTorch provides when processing batches (in our case the approximation of multiple evaluations from ${\phi }_n\left({\omega }_i\right)$ with NNs). However, the amount of evaluations of ${\phi }_n\left({\omega }_i\right)$ that Algorithm A2 can process at once is limited by the number of available processors. Since the calculation time for the evaluation of an optimisation step ${\phi }_{n+1}\left({\omega }_i\right)$ increases with finer meshes, the evaluation of the approximation of all gradient steps ${\phi }_{n+1}\left({\omega }_1\right),\dots ,{\phi }_{n+1}\left({\omega }_S\right)$ at once results in a processing time advantage for the NN. It is to be expected that this time saving increases with the number of examples S.

As mentioned at the beginning of the experiment, we expect to achieve good results close to the mean of $g={(0,-5000)}^T$. In order to obtain a more general view of the quality of Algorithm 1, we have compiled a selection of extreme cases for the distribution of g (e.g., evaluations from g that deviate strongly from ${(0,-5000)}^T$) in Figure 8 and

Table 3. Comparison of metrics between Algorithms A1 and 1.

Method	Applied Load g	Converges after ${\mathit{n}}_{\mathbf{final}}$	Evaluation of ${\mathit{J}}^{\mathit{\epsilon }}\left({\mathit{\phi }}_{{\mathit{n}}_{\mathbf{final}}}\right)$	Compliance ${\int }_{{\mathbf{\Gamma }}_{\mathit{g}}}\mathit{g}·\mathit{u}\left(\mathit{\phi }\right)\mathbf{d}\mathit{s}$
Algorithm A1	${(0,-5000)}^T$	538	1475.39	1173.22
Algorithm 1	${(0,-5000)}^T$	441	1206.49	1148.66
Algorithm A1 (see Figure 8a)	${(2632.16,-4251.09)}^T$	473	1487.39	1416.46
Algorithm 1 (see Figure 8a)	${(2632.16,-4251.09)}^T$	517	1475.05	1405.39
Algorithm A1 (see Figure 8b)	${(-733.23,-4945.95)}^T$	457	1115.66	1062.41
Algorithm 1 (see Figure 8b)	${(-733.23,-4945.95)}^T$	497	1049.91	1042.58
Algorithm A1 (see Figure 8c)	${(-1099.92,-4877.25)}^T$	470	1042.18	992.28
Algorithm 1 (see Figure 8c)	${(-1099.92,-4877.25)}^T$	447	1048.44	998.34

. The figure shows the sequence of ${\phi }_n$ in hundreds of steps as well as the final distribution of material (${\phi }_{100},{\phi }_{200},\dots ,{\phi }_{n_{\mathrm{final}}}$) for the specific loads g. Table 3 indicates a noticeable saving in calculation time, but there is no guaranteed improvement in the results. In particular, when the topology “collapses” (i.e., the NN cannot generalise to the input data with strong deviations from the training data), the application of the CNN leads to worse results. Nevertheless, it can be seen that the NN extension gives the algorithm a greater robustness against porous fragments (see Figure 8b) in the optimisation of the topology and thus a higher stability against collapsing of the topology in the optimisation can be assumed. Finally, a critical aspect to be mentioned is the step size $c_m$. The time at which Equation (22) is applied and which is controlled by $c_m\in \mathbb{N}$ has a crucial impact on the viability or “compatibility” between the state of the optimisation procedure and the CNN. In some cases, a $c_m$ that is too small or too large can lead to the collapse of the topology, i.e., the topology deteriorates and does not recover. For the reliable use of Algorithm 1, a method for controlling $c_m$ would have to be devised.

4.2. TLSTM Examples

As in Section 4.1.3, randomly generated training data should be used in the following experiment with the derived LSTM transformation of Equation (31). The data tuples consist of input and output from ${\mathcal{N}}_{\mathrm{LSTM}}$ according to Example 4.

4.2.1. Sampling the Data

Again, we assume an expected load $g={(0,-5000)}^T$ and a random rotation characterised by a truncated normal distribution with bounds $[-\frac{\pi }{2},\frac{\pi }{2}]$, standard deviation $0.3$ and mean 0. Using the optimiser in Algorithm A1, $S\in \mathbb{N}$ sample paths of gradient steps $\phi \left(g\right)$ and solutions of the state equation $u\left(g\right)$ are generated by drawing S realisations of g. In contrast to Section 4.1.1, this time we do not only store every $T\in \mathbb{N}$ iteration step of ${\phi }_n$ and $u_n$ but instead store all iteration steps of the optimisation of Algorithms A1. Afterwards, these are merged into disjoint subsets, each consisting of a sequence of T iteration steps. Thus, the feature or the sequence ${\phi }_{n+1},\dots ,{\phi }_{n+T}$ is the label for the assembled sequence from ${\phi }_{n-T},\dots ,{\phi }_n$ and $u_{n-T+1},\dots ,u_{n+1}$. More precisely, with

$$\left(\left[\begin{array}{ccc}{\phi }_{n-T}& u_{n-T+1}^x& u_{n-T+1}^y\\ & ⋮& \\ {\phi }_n& u_{n+1}^x& u_{n+1}^y\end{array}\right],\left[\begin{array}{c}{\phi }_{n+1}\\ ⋮\\ {\phi }_{n+T}\end{array}\right]\right),$$

for $0\le n\le N-(T-1)$, a total of $S(\lfloor \frac{N}{k}\rfloor -1)$ tuples are stored in an unsorted dataset ${\mathcal{D}}_{\mathrm{LSTM}}$.

4.3. TLSTM Predictions

After training the TLSTM from Equation (29) with the data set ${\mathcal{D}}_{\mathrm{LSTM}}$ generated in Section 4.2.1, we wish to investigate its predictive ability of ${\mathcal{N}}_{\mathrm{LSTM}}^T\left({\phi }_n\right):={\mathcal{N}}_{\mathrm{LSTM}}^T({\phi }_{n-T},\dots$$,{\phi }_n,u_{n-T+1},\dots ,u_{n+1};\theta )$ compared to the real sequence ${\left({\phi }_{n+T}\right)}_n$. For this purpose we have visualised both sub-sequences for $T=10$ in Figure 9 using the iteration sequence generated for a load $g=(0,-5000)$. The distorted topology in the first forecasts is striking. This can be attributed to the comparatively low weighting of training data in which the distribution of the material is constant ${\phi }_n\left(x\right)=0,5$ for all $x\in D$ or almost constant. This forces us to choose a correspondingly high $c_m\in \mathbb{N}$ in Algorithm 3. Furthermore, it can be seen that especially in early phases of the partial sequence in which the change $∥{\phi }_n-{\phi }_{n+1}∥$ is very high, ${\mathcal{N}}_{\mathrm{LSTM}}^{10}\left({\phi }_n\right)$ provides a better forecast from a visual perspective, i.e., the topology ${\mathcal{N}}_{\mathrm{LSTM}}^{10}\left({\phi }_n\right)$ has already converged further than the target image ${\phi }_{n+10}$. Since the topologies on the finer meshes no longer show any major visual changes and therefore the differences in the predictions are no longer recognisable, we have decided not to present them at this point.

The architecture of the TLSTM allows the length of the input sequence as well as the output sequence to be chosen independently of the training data. The expected consequence is a decrease in prediction quality. Despite this, Figure 10 depicts the prediction results of Equation (29) with unchanged input sequence ($T_{en}=5$) and output sequence of length 40. Since a shorter output sequence ($T_{des}=10$) is used to train Equation (29), the results of the longer output sequence indicate that $N_{\mathrm{LSTM}}$ has indeed learned to predict a gradient step for the given setting and that the training data from Section 4.2.1 describe the problem correctly.

Analogous to Section 3.1.1, the intention behind the construction of ${\mathcal{N}}_{\mathrm{LSTM}}$ is to replace the gradient step in reference Algorithm A1 with Example 3. Algorithm 3 describes the integration of the LSTM prediction. When evaluating ${\mathcal{N}}_{\mathrm{LSTM}}^T$, only the last iteration step ${\phi }_{n+T}$ of the predicted sequence is projected back to the current mesh ${\mathcal{T}}_m$ in order to save computational resources. We examine the performance of Algorithm 3 in the following deterministic experiment. As for the TCNN above, the beneficial performance of a single-gradient prediction transfers to the stochastic setting since it consists of a Monte Carlo estimator with $N\in \mathbb{N}$ samples in each step. It is hence not necessary to examine this in more detail.

4.4. Deterministic Bridge Optimisation

The motivation for the design of the TLSTM architecture was that the information contained in the time series of ${\phi }_n$ and $u_n$ could ideally lead to an improvement in the forecast capabilities of the NN. This can be investigated as in Section 4.1.4 by calculating the optimal topology for different loading scenarios for g by Algorithm 3. Again, all results are based on the ten-fold averaged performance of the algorithms for each load g. Figure 11 shows the results in the setting similar to Appendix A.1 and compares the respective metrics. Analogous to Section 4.1.3 the sequence ${\phi }_n$ of the optimisation by Algorithm 3 is also more resilient to porous fragments in the structures than the reference optimisation procedure. In general, the pictures of ${\phi }_n$ hardly differ between Algorithms 1 and 3. Hence, apart from Figure 11, no further visualisations are presented. It should also be noted that the optimisation by Algorithm 3 is much less stable than the optimisation using ${\mathcal{N}}_{\mathrm{CNN}}$ from Equation (23). This becomes apparent when the structure collapses which was the case in each of our test runs if the $c_m\in \mathbb{N}$ chosen was too small in Algorithm 3. Furthermore, it could be observed that the convergence criterion of Equation (A1) was not reached after applying ${\mathcal{N}}_{\mathrm{LSTM}}$ because ${\phi }_n$ diverged too far from the actual minimum on ${\mathcal{T}}_m$. In conclusion, the stability of Algorithm 3 is even more dependent on $c_m$ than it is with Algorithm 1, which renders parameter calibration more difficult. However, the metrics in Figure 11 show that Algorithm 3 converges faster than Algorithm A1 and often achieves better results.

One conspicuous feature is the high fluctuation in the calculation time per iteration in the applications of ${\mathcal{N}}_{\mathrm{LSTM}}(·;\theta )$ during the optimisation when compared to Algorithm 1. On the one hand, this is due to the comparatively high complexity of the TLSTM. In this context, high complexity means a high dimension $d\in \mathbb{N}$ of the parameter vector $\theta \in {\mathbb{R}}^d$. On the other hand, the main driver of the higher calculation time is the transformation given by Equation (31) since on an algorithmic level the entire input sequence ${\phi }_n,{\phi }_{n+1},{\phi }_{n+2},{\phi }_{n+3},{\phi }_{n+4}$ has to be stored and transformed. In general, it becomes apparent at this point that the transformations in Equations (21) and (31) are the critical aspects that compromise the performance of Algorithms 1 and 3.

Table 4. Comparison of metrics between Algorithms A1, 1 and 3.

Method	Applied Load g	Computation Time	Converges after ${\mathit{n}}_{\mathbf{final}}$	Evaluation of ${\mathit{J}}^{\mathit{\epsilon }}\left({\mathit{\phi }}_{{\mathit{n}}_{\mathbf{final}}}\right)$	Compliance ${\int }_{{\mathbf{\Gamma }}_{\mathit{g}}}\mathit{g}·\mathit{u}\left(\mathit{\phi }\right)\mathbf{d}\mathit{s}$
Algorithm A1	${(0,-5000)}^T$	4 min 43 s	538	1475.39	1173.22
Algorithm 1	${(0,-5000)}^T$	4 min 05 s	441	1206.49	1148.66
Algorithm 3	${(0,-5000)}^T$	4 min 34 s	425	1209.24	1151.27
Algorithm A1	${(2632.16,-4251.09)}^T$	4 min 31 s	473	1487.39	1416.46
Algorithm 1	${(2632.16,-4251.09)}^T$	4 min 14 s	517	1475.05	1405.39
Algorithm 3	${(2632.16,-4251.09)}^T$	4min 34 s	436	1209.24	1151.27
Algorithm A1	${(-733.23,-4945.95)}^T$	4 min 50 s	457	1115.66	1062.41
Algorithm 1	${(-733.23,-4945.95)}^T$	4 min 21 s	497	1049.91	1042.58
Algorithm 3	${(-733.23,-4945.95)}^T$	4 min 41 s	484	1082.62	1031.62
Algorithm A1	${(-1099.92,-4877.25)}^T$	5 min 11 s	470	1042.18	992.28
Algorithm 1	${(-1099.92,-4877.25)}^T$	4 min 42 s	447	1048.44	998.34
Algorithm 3	${(-1099.92,-4877.25)}^T$	4 min 43 s	542	1041.28	992.16
Algorithm A1	${(-3197.16,-3844.24)}^T$	2 min 17 s	365	845.93	805.94
Algorithm 1	${(-3197.16,-3844.24)}^T$	2 min 13 s	346	852.36	811.15
Algorithm 3	${(-3197.16,-3844.24)}^T$	2 min 59 s	366	848.52	808.40

compares the performance of the three presented algorithms. In overall terms, Algorithm 3 achieves better compliance, whereas Algorithm 1 stands out owing to its shorter calculation time.

5. Discussion and Conclusions

The objective of this work is to devise neural network architectures that can be used for efficient topology optimisation problems. These tasks are computationally burdensome and typically are inevitably carried out with a large number of optimisation steps, each requiring (depending on the chosen method) the solution of state and adjoint equations to determine the gradient direction. Instead of learning a surrogate for state and adjoint equations, we present NN architectures that directly predict this gradient, leading to very efficient optimisation schemes. A noteworthy aspect of our investigation is the consideration of uncertainties of model data in a risk-averse optimisation formulation. This is a generalisation of the notion of “loading scenarios” that are commonly used in practice for a fixed set of parameter realisations. With our continuous presentation of uncertainties in the material and of the load acting on the considered structure, the robustness of the computed design with respect to these uncertainties can be controlled by the parameter of the CVaR used in the cost functional. Since computations with uncertainties require a substantial computational effort, our central goal is to extend the algorithms used in [1,2] by introducing appropriate NN predictions, reducing the iteration steps required. In contrast to other machine learning approaches, our aim is to achieve this even for adaptively adjusted finite element meshes since this has proven to be crucial for good performance in previous work. For this to function, an underlying sufficiently fine reference mesh is assumed for the training data and the prediction. Moreover, in contrast to other NN approaches for this problem, we consider the evolution of a continuous (functional) representation of a phase field determining the material distribution.

Ideally, the NN architectures should hasten the deterministic topology optimisation problem and consequently the risk-averse optimisation under uncertainties. This is achieved in Section 3.1 by embedding a CNN in the optimisation for both the deterministic and the stochastic setting.

The observed numerical results for a common 2D bridge benchmark are on par with the reference method presented in [1]. However, the gradient step predicted by the NN architectures allows for significantly larger iteration steps, rendering the optimisation procedure more efficient. This directly transfers to the Monte-Carlo-based risk-averse optimisation under uncertainties as defined by Algorithm 2 since the samples for the statistical estimation are obtained with minimal cost. In addition to the CNN, a second architecture is illustrated in terms of an LSTM. This generally leads to a better quality of the optimisation and is motivated by the idea that a memory of previous gradients may lead to a more accurate prediction of the next gradient step. However, it comes at the cost of longer computation times due to the transformation between the different adapted computation meshes (see Section 4.4). Hence, a substantial performance improvement could be achieved by reducing the complexity of the transformations of Equations (21) and (31).

There are several interesting research directions from the presented approach and observed numerical results. Regarding the chosen architectures, an interesting extension would be to consider graph neural networks (GNN) since there, the underlying mesh structure is mapped directly to the NN. Consequently, the costly transfer operators from current mesh to reference mesh of the design space could be alleviated, removing perhaps the largest computational burden of our approach. Moreover, transformer architectures have probably superseded LSTMs and it would be worth examining this modern architecture in the context of this work.

The loss function used in the training also leaves room for improvements. For example, instead of the simple mean squared error used here, one could approximate the objective functional of Equation (5) directly in the loss function. Regarding the training process, there are modern techniques to improve the efficiency and alleviate over-fitting such as early stopping, gradient clipping, adaptive learning rates and data augmentation as discussed in [25]. Moreover, transfer learning in a limited-data setting could substantially reduce the amount of training data required.

This work mainly serves as a proof of concept for treating the considered type of optimisation problems with modern NN architectures. An important step towards practicability is the further generalisation of this model, e.g., to arbitrary problems (with parameterised boundary data and constraints), determined by descriptive parameters drawn from arbitrary distributions according to the problem at hand. Moreover, the models presented here can be used as a basis for theoretical proofs (e.g., regarding the complexity of the representation) and further systematic experiments.