Improving the Accuracy of the Deep Energy Method

1. Introduction

With the advancements in computational technology and complexity in design, engineers seek to predict the component’s performance before it can be manufactured. These predictions involve complex boundary value problems (BVP), which require solutions to partial differential equations (PDEs). Typically, these PDEs are solved through numerical approximations. In solid mechanics, conventional methods to solve PDEs include mesh-free methods^1,2, finite element analysis (FEA)³, and isogeometric analysis^4,5. However, these methods are poised with challenges, such as obtaining robust and accurate solutions for ill-posed, high-dimensional, or coupled problems, to name a few.

Machine learning methods are rapidly developing as an alternative to traditional approaches to minimize or, in some cases, eliminate the issues mentioned above. Liu et al.⁶ have discussed FEA’s evolution, future, and transition to machine learning applications in detail. In the case of solid and fluid mechanics, machine learning models can broadly be classified into data-driven models^{6,7,16,17,8–15} and physics-informed neural networks (PINNs)^18–25. In data-driven models, data from experimental and computational results are used to train the models. Although data-driven models can capture complex physical phenomena, the amount of data required to train the models, data quality, and data generation process impede its widespread implementation²⁶. Another application of data-driven models in mechanics is to provide an automated framework for predicting constitutive models for material behavior, as discussed by Flaschel et al.²⁷, Yang et al.²⁸, and Chen²⁹. The data-driven models can also take images to identify the material distribution and defect characterization³⁰. To reduce dependency on a large data set, Saha et al.³¹ developed a model (hierarchical deep learning neural network (HiDeNN)) with universal approximation capability to interpolate between data points for extremely nonlinear relationships.

In contrast to data-driven models, PINNs use a deep neural network (DNN) to train the model based on physical laws at a set of points in the domain. PINNs with automatic differentiation were first introduced by Raisi et al.¹⁸. They have several advantages over data-driven models as they do not require data labeling and require no or minimal preliminary datasets compared to purely data-driven models. Two types of PINN models are rapidly developing in solid mechanics: deep collocation method (DCM)^23–25,32 and deep energy method (DEM)^33–35. In DCM, residuals of strong form define the loss function at points sampled from the physical domain, corresponding boundary, and initial conditions (collocation points). This method has been extended and improved by researchers^36,37. Haghighat et al.³⁸ used a PINN technique to solve the coupled flow and mechanical deformation equations in porous media for single-phase and multiphase flows. Improving the application of PINNs to irregular geometries, Gao et al.³⁹ developed PhyGeoNet, a convolutional neural network model (CNN) that could learn solutions for parametric PDEs on irregular domains using elliptic coordinate mapping. Chen et al.⁴⁰ generalized the PINNs to discover governing PDEs from sparse and noisy data. Though the DCM has been widely successful in predicting outcomes, the approach requires higher-order gradients, which can be computationally expensive.

In the case of DEM, the system’s potential energy (PE), expressed as a loss function, is minimized to predict the system’s displacements. Compared to classical PINNs, DEM has advantages as it relies only on first-order differentiation to train the neural network and on accurate numerical integration. The DEM method can be extended to the problems where an energy function exists and reduces dependencies on PDEs of the base function. E and Yu⁴¹ first proposed the idea of DEM using the Ritz method to solve variational problems. Nguyen-Thanh et al.³⁴ expanded the work presented by E and Yu to 2D and 3D material models. Implementation of DEM for several BVPs has been reported in the literature^33,35,42. Extending this work, Fuhg and Bouklas⁴³ demonstrated that DEM and DCM fail to accurately resolve displacement and stresses at stress concentration regions. They proposed modified DEM (mDEM) to overcome this issue and added stresses to the loss function. However, it is worth noting that including stresses in the loss function requires a second-order differential equation. The modification allowed the successful resolution of stresses around the stress concentration regions. Abueidda et al.⁴⁴ enhanced the model by developing PINN that combines residuals of the strong form and the system’s potential energy. The proposed formulation yielded a loss function with multiple loss terms. As a result, the coefficient of variation weighting scheme was also introduced in the loss function to assign the weight dynamically and adaptively for each term. He et al. extended the applicability of DEM to plasticity⁴⁵ and graph convolutional network⁴⁶. Similar to DEM, Sun et al.⁴⁷ used DNN to develop surrogate models for fluid flows using a loss function that incorporated initial boundary conditions and governing PDEs. One of the disadvantages of PINNs over classical FEA is the time required to train the neural network model. However, users can use transfer learning^48,49, where instead of retraining the model from scratch, the network can either be partially retrained, trained using prior weights and biases, or use a trained model to predict mechanics for similar geometry.

Though DEM has been successfully used to predict deformations in bending, similar accuracies were not obtained for compression and tension loadings, even for a simple 2D rectangular linear elastic plate. Figure 1 shows the deformation under compression loading in the y-direction using the architecture proposed by Nguyen-Thanh et al.³⁴. The geometry and boundary conditions for the plate are shown in Figure 7. From Figure 1, we can notice that even though symmetric pressure is applied on the plate, the resulting displacements obtained from DEM are not symmetric. Thus, there is a need to improve the accuracy of DEM under generalized loading conditions.

In this paper, two modifications are proposed to improve the accuracy of DEM while retaining dependency only on first-order derivatives. Firstly, random Fourier feature mapping is introduced in the multilayer perceptron (MLP) model to reduce bias that may be present during the model’s training. Secondly, a two-loop architecture is proposed to obtain appropriate hyperparameters for a given geometry. Recent studies have demonstrated that optimizing hyperparameters in large and multi-layered models impedes scientific progress⁵⁰. This process can be challenging for DEM due to the high interaction between independent hyperparameters, large search space, difficulty in identifying objective functions, and non-convex relationships with the objective functions. Thus, the effect of different hyperparameters in loss function and the effect of loss function on displacement are also studied to find a generalized search space.

Tuning these hyperparameters through manual trial and error is highly time-consuming. As a result, researchers have employed grid search, random search, and optimization algorithms to obtain the best combination of hyperparameters^22,50–53. Still, studies on optimizing hyperparameters and architectures of PINNs for solving mechanics problems remain uncommon. Wang et al.⁵⁴ proposed Auto-PINN, a framework for employing Neural Architecture Search (NAS) techniques for PINNs to solve heat transfer problem. However, unlike Wang et al.⁵⁴, we found that not all hyperparameters can be decoupled for DEM.

In addition, the L₂ norm is commonly employed as an objective function to obtain optimized hyperparameters. This approach requires finding a solution to the BVP beforehand. Since hyperparameters can be problem-specific, using the L₂ norm as an objective function increases the dependency of PINNs on existing solutions, reducing their generalizability. The current research demonstrates that the principle of minimum potential energy can be used as an alternative objective function to obtain optimized hyperparameters. This approach mitigates the need for prior knowledge of the solution. The examples presented in the current research are based on 2D linear elastic plane stress problems. Even though the examples solved in this study involve 2D linear elastic material models, the proposed approach can be extended to a variety of DEM problems (for example, to 3D and nonlinear material models). Random Fourier feature (RFF) mapping (discussed in Section 3) is introduced to improve the accuracy of DEM further. RFF mapping helps in reducing bias during learning towards high frequencies. The effect of reduction in bias becomes prominent in complex geometries and loading conditions.

The paper is divided into eight sections. Section 2 and Section 3 describe the formulation of DEM and modifications made to DEM (m-DEM). Preliminary observations relating the accuracy of displacements obtained from DEM to the loss function and hyperparameters are discussed in Section 4. Based on initial observations, a search space is defined and followed by hyperparameter optimization for the elastic rectangular plate in Section 5. The transferability of optimized hyperparameters obtained in Section 5 to different geometry and load cases is discussed in Section 6. The effect of adding RFF mapping in MLP is shown in Section 7. Finally, conclusions are presented in Section 8.

2. DEM formulation

The DEM consists of a multilayer perceptron (MLP) neural network, as shown in Figure 2. In general, MLP consists of sets of neurons arranged in layers. Each neuron is connected to all the neurons in adjacent layers. The number of layers (N) and neurons in each layer are predefined according to specific applications. In DEM, the number of neurons in the first and last layer (also known as input and output layers, respectively) is determined by the system’s dimension (two in the case of 2D and three in 3D). The programmer predefines the number of neurons in the layers between the input and output layers (also called hidden layers).

Mathematically, the value of each neuron connecting the previous layer ($n_{l-1})$ is calculated using a linear combination of weights (w) and biases (b), as shown in Equation 1. An activation function (φ: ℝ→ ℝ) then operates on the calculated value to determine the value passed by the neuron.

$$\widehat{{\mathit{z}}^l}=\phi \left({\mathit{w}}^l{\mathit{z}}^{l-1}+{\mathit{b}}^l\right)$$

(1)

The values of weights and biases are determined such that the loss function ($\mathcal{L}\left(\widehat{z^l}(w,b)\right)$ or $\mathcal{L}\left(\mathsf{\Theta }\right)$) is minimized. In a so-called backpropagation procedure, an optimization algorithm can be used to iterate over the loss function to predict values of weights and biases such that:

$${\mathsf{\Theta }}^{*}={arg}_{\mathsf{\Theta }}min\mathcal{L}\left(\mathsf{\Theta }\right)$$

(2)

Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm is used in the examples shown in the following sections to minimize the loss function and demonstrate the performance of DEM. However, other algorithms can also be used. The backpropagation capabilities of PyTorch are used to calculate gradients of the loss function.

In the above description, parameters like the number of layers and of neurons in each layer need to be defined before MLP can be trained. These parameters, called hyperparameters, define the architecture of the MLP and control the learning process of the model. It is essential to obtain optimized values for the hyperparameters to achieve a minimum value for the loss function. The number of hyperparameters and their values vary according to the problem statement. Six hyperparameters based on the architecture of DEM were determined. These are the total number of layers, number of neurons in hidden layers, activation function, standard deviation for predicting weights and biases, learning rate for the L-BFGS algorithm, and the number of epochs for the L-BFGS algorithm.

The loss function for DEM is defined as the system’s potential energy (Equation 3), where U is the system’s internal energy, W is the energy due to externally applied load, $\mathit{v}$ is a kinematically admissible displacement field, H is the space of admissible functions, $\psi$ is strain energy density function, ${\mathit{f}}_b$ are internal forces of the body B, and $\tilde{\mathit{t}}$ are traction forces applied to the area $\text{∂}{B}_t$ (Figure 3).

$$\mathcal{L}\left(\mathit{v}\right)=\underset{\mathit{v}\in H}{\min }\mathsf{\Pi }\left(\mathit{v}\right)=\underset{\mathit{v}\in H}{\mathrm{m}\mathrm{i}\mathrm{n}}(U-W)=\underset{B}{\int }\psi \left(\mathit{ϵ}\right)dV-\underset{B}{\int }{\mathit{f}}_b.\mathit{v}dV-\underset{\text{∂}{B}_t}{\int }\tilde{\mathit{t}}.\mathit{v}dA$$

(3)

The selection of loss function is based on the principle of minimum potential energy: out of all the possible kinematically admissible deformations, a conservative structural system assumes displacements that minimize total potential energy at equilibrium.

3. Modifications to DEM

3.1. Random Fourier feature mapping

Fourier transformation converts functions dependent on space or time domain into functions dependent on spatial or temporal frequency. This transformation is widely used in mathematics and digital image processing to simplify calculations. Tancik et al.⁵⁵ demonstrated that the standard MLP fails to learn high frequencies, and Fourier feature mapping can be used to mitigate this bias. Wang et al.⁵⁶ expanded the work of Tancik et al.⁵⁵ to PINNs, discussed the limitations of PINNs via the neural tangent kernel theory (NTK), and illustrated how PINNs are biased towards learning along the dominant eigenvectors of their limiting NTK. Correspondingly, they have demonstrated that using Fourier feature mappings with PINNs helps modulate the frequency of the NTK eigenvectors.

Figure 4 shows the implementation of the Fourier feature mapping. The input layer (coordinates of the training points from the physical domain) is passed through a Fourier mapping γ before passing through MLP. The function γ maps the input coordinates to a higher-dimensional hypersphere using a set of sinusoids:

$$\gamma =[a_1\cos \left(2\pi h_l^{T}X\right),a_1\sin \left(2\pi h_l^{T}X\right),\dots \dots ,a_m\cos \left(2\pi h_m^{T}X\right),a_m\sin \left(2\pi h_m^{T}X\right)]$$

(4)

where h_i’s are the Fourier basis frequencies used for approximation and a_i represent the corresponding coefficients in the Fourier series. Here, we adopted Gaussian random Fourier feature (RFF) mapping. Each entry in h is sampled from a normal distribution N (0, σ²), where σ² is the variance. σ is a hyperparameter that needs to be specified for each problem. Adding RFF mapping helped improve the solution’s accuracy even with optimized hyperparameters. The details on the effect of RFF mapping are discussed in Section 7.

3.2. Tuning of hyperparameters

Like most MLPs, the accuracy of DEM solutions depends on selecting appropriate hyperparameters. We will see in Section 4 that the hyperparameters are highly interdependent, and the search space is highly non-convex. In addition, even for a simplified linear elastic material model, getting stuck in local minima during hyperparameter tuning and obtaining incorrect solutions is relatively easy. A two-loop architecture (as shown in Figure 5) is proposed to overcome these issues and simplify the hyperparameter tuning. The two loops are referred to as the inner loop (solves DEM with Fourier feature mapping) and the outer loop (used for hyperparameter tuning).

Three essential things must be considered to initialize a successful hyperparameter optimization: search space, optimization algorithm, and objective function. The selection of search space is discussed in detail in Section 4. Since the search space is non-convex and consists of continuous and discrete variables, a Bayesian optimization algorithm in the HyperOpt library (tree pazen estimator: TPE) was selected for the current study⁵⁷. It is worth noting that the algorithms used for optimizing hyperparameters are an evolving field, and new algorithms are being rapidly developed. However, the selection of the best optimization algorithm for DEM is outside the scope of the current study.

Minimization of PE (Equation 3) was chosen as the objective function for hyperparameter tuning. The selection was based on the principle of minimum potential energy. The advantages of choosing minimal PE are minimized additional calculations required during the tuning process and no required prior knowledge of the solution. Thus, the framework can be extended to any problem statement that can be solved using DEM, and hyperparameters can be turned to specific BVP if required without knowledge of the exact solution. The algorithm for the tuning process is described in Figure 6. A BVP is defined along with search space and profiles for hyperparameters. The hyperparameters are initialized randomly based on the search space and supplied to DEM. DEM then solves for weights and biases to obtain displacements for minimum PE. Once convergence has been reached, the value for minimum PE is supplied to the outer loop, where an observation set D is formed. A surrogate model is formed using D to predict the next set of hyperparameters. The new hyperparameters are supplied to the inner loop, which runs DEM to solve for BVP and supplies minimum PE. D is updated based on new results. The cycle continues until the termination condition is satisfied. We used two termination conditions for examples discussed in Section 5. The first condition asks the algorithm to be terminated if no improvement in minimum PE is observed for thirty consecutive iterations. This condition also specifies the minimum number of iterations that need to be undertaken by the algorithm. The second condition asks for the algorithm to be terminated after fifty iterations. The condition can be varied according to the maximum number of iterations that should be run for hyperparameter tuning.

The initial TPE algorithm was proposed by Bergstra et al.⁵⁸. This approach employs both estimator models p(x|y) and p(y) to reduce computations. It transforms the prior distribution of the search space into truncated Gaussian distribution and modifies the posterior distribution based on observations. It then sorts the target value (y) and divides it into two, as shown in Equation 5.

$$p(x|y)=\left\{\begin{array}{c}l\left(x\right)\\ g\left(x\right)\end{array}\right.\begin{array}{c},y<y^{*}\\ ,y\ge y^{*}\end{array}$$

(5)

where $y^{*}$ is the boundary used to segregate the target value, $l\left(x\right)$ is the density formed by observations that are less than $y^{*}$, and g(x) is the density formed by the remaining observations. The TPE algorithm chooses $y^{*}$ to be some quantile ζ of the observed y values so that p(y < y^∗ ) = ζ. The next promising point is chosen based on the Expected Improvement (EI) criterion⁵⁹.

$${EI}_{y^{*}}\left(x\right)={\int }_{-\infty }^{y^{*}}\left(y^{*}-y\right)p(y|x)dy={\int }_{-\infty }^{y^{*}}\left(y^{*}-y\right)\frac{p(y|x\left)p\right(y)}{p\left(x\right)}dy$$

(6)

4. Pre-experiments and observations

The dependency of PE on hyperparameters can be best understood using a simplified BVP. Thus, an elastic plate subjected to compression and bending loads (shown in Figure 7) is shown as an example. The plate is 4 mm in length and 1 mm wide. It is composed of an elastic material having young’s modulus of 1000 MPa and Poisson’s ratio of 0.3. The plate is subjected to uniform compression and bending loads, as shown in Figure 7b,c. For compressive load, a uniform compressive pressure of 50 MPa was applied to the right, while the left end was fixed in all degrees of freedom. In the case of bending loads, a uniform transverse shear load of 10 MPa^, forcing the plate to move downwards, was applied on the right, and the left end was fixed in all degrees of freedom. The strain energy density of the plate is given by:

$$\psi \left(\mathit{ϵ}\right)=\mu \left(\mathit{ϵ}:\mathit{ϵ}\right)+\frac{1}{2}\lambda {\left(\mathit{ϵ}:1\right)}^2$$

(7)

where $\mu$ and $\lambda$ are Lame constants. $\mathit{ϵ}$ is the strain tensor for small deformations, which is related to displacements $\mathit{u}$ by:

$$\mathit{ϵ}=\frac{1}{2}\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)$$

(8)

The relationships between the Lame constants, young’s modulus (E) and Poisson’s ratio (ν) are given by equation 9.

$$\mu =\frac{E}{2(1+v)}$$

(9a)

$$\lambda =\frac{Ev}{(1+v)(1-2v)}$$

(9b)

Differentiation of strain energy density with respect to strain tensor provides Cauchy stress tensor $\mathit{\sigma }$, which can be differentiated again to find the stiffness tensor C.

$$\mathit{C}=\frac{\partial \mathit{\sigma }}{\partial \mathit{ϵ}}=\frac{{\partial }^2\psi \left(\mathit{ϵ}\right)}{\partial \mathit{ϵ}\partial \mathit{ϵ}}$$

(10)

For a 2D plane stress condition the relation between Cauchy stress tensor and strain tensor can be reduced to:

$$\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}2\mu +\lambda & \lambda & \lambda \\ \lambda & 2\mu +\lambda & \lambda \\ 0& 0& 2\mu \end{array}\right]\left\{\begin{array}{c}{ϵ}_{xx}\\ {ϵ}_{yy}\\ {ϵ}_{xy}\end{array}\right\}$$

(11)

4.1. Interactions effect

Sensitivity analysis can be used to determine the impact of independent inputs (hyperparameters) on outputs (PE under compression and bending). A five-level Taguchi design of experiment (DOE) and two-level half-factorial DOE were employed to obtain the dataset for sensitivity analysis. The values of hyperparameters used to generate a five-level Taguchi DOE are shown in Table 1. Levels 1 and 5 from Table 2 were utilized to design the two-level half-factorial DOE.

After running m-DEM, the minimum potential energy corresponding to each set of hyperparameters was collected, and the data were analyzed using Minitab. Similar results were obtained for compression and bending loads. The results obtained from compression loading conditions are discussed in this section.

The analysis of the data was conducted using Minitab. Figure 8 and Figure 9 show the main effect and interaction plots between hyperparameters. Figure 8 illustrates that the mean response (minimum PE) varies for each hyperparameter. We also notice that the minimum PE obtained from DEM does not follow the same trends across hyperparameters, and it varies least with the variation in activation function.

Table 2. Values used to generate data for sensitivity analysis.

Variable	Level 1	Level 2	Level 3	Level 4	Level 5
Learning Rate	0.1	0.368	1	2.718	7.389
Neurons	20	52	70	96	120
Standard Dev (DNN)	0.0001	0.001	0.01	0.1	1
Standard Dev (RFF)	0.0001	0.001	0.01	0.1	1
Total number of layers	3	4	5	6	7
Activation function	rrelu	relu	celu	sigmoid	tanh

Table 2. Values used to generate data for sensitivity analysis.

Variable	Level 1	Level 2	Level 3	Level 4	Level 5
Learning Rate	0.1	0.368	1	2.718	7.389
Neurons	20	52	70	96	120
Standard Dev (DNN)	0.0001	0.001	0.01	0.1	1
Standard Dev (RFF)	0.0001	0.001	0.01	0.1	1
Total number of layers	3	4	5	6	7
Activation function	rrelu	relu	celu	sigmoid	tanh

Figure 9 helps us visualize significant interaction effects between the hyperparameters even though each can be varied independently. Due to the high interaction, multiple hyperparameters should be varied simultaneously for tuning the hyperparameters. This process can be very tedious and time-consuming if done manually. Alternatively, computational algorithms can be used to speed up the process.

In addition to plotting the main effects plot and interaction plot, the normality of the loss function obtained from two DOE was tested using the Anderson-Darling test. Figure 10 shows the probability plot for the loss function. Based on the test results, we can notice that the data is not normally distributed (p-value <0.005). We can also see the high probability of obtaining a loss function near -4.5. The importance of this observation is discussed later in this section.

Loss functions from five hundred DEM iterations were collected to further study the six hyperparameters’ effect. The hyperparameters were varied using the random algorithm present in HyperOpt. The minimum potential energy obtained from the iterations plotted against hyperparameters is shown in Figure 11. However, from this figure, we do not see any trend between hyperparameters and the loss function.

4.2. Sensitivity of loss functions and displacements

Figure 11 shows that the loss function for most iterations lies in a narrow range of -5 to -4.5. Displacements in the y-direction for two randomly picked iterations are shown in Figure 12. From this figure, we can see that even though the variation in loss function is relatively small, neither provides the correct displacement prediction. Previously we had observed that the loss function is not normally distributed, and without a trend, it is not easy to tune hyperparameters manually. This observation further strengthens the requirement for a systematic approach to tuning hyperparameters.

5. Tuning hyperparameters for different load cases

In the previous section, we concluded that all six hyperparameters are interdependent. However, from the main effects plot (Figure 8), we can estimate that the activation function has a lesser impact on the loss function than other parameters. From random five hundred iterations, we see that the median loss function for rrelu and sigmoid is less than other employed activation functions. As a result, the activation function could be fixed to one of them, and here we select rrelu. Further, the number of layers was fixed to five. The other four hyperparameters varied according to the sensitivity profiles and ranges provided in Table 3.

Table 3. Parameter ranges and profiles.

Hyperparameter	Sensitivity profile	Range
Learning Rate	loguniform	Exp(0-2)
Neurons	quniform	20-120
Standard Dev (DNN)	uniform	0-1
Standard Dev (RFF)	uniform	0-1

Table 3. Parameter ranges and profiles.

Hyperparameter	Sensitivity profile	Range
Learning Rate	loguniform	Exp(0-2)
Neurons	quniform	20-120
Standard Dev (DNN)	uniform	0-1
Standard Dev (RFF)	uniform	0-1

Sensitivity profile and ranges together define the search domain. While the range provides the upper and lower bounds for each hyperparameter, the sensitivity profile influences how values are selected within the range. A profile of loguniform returns values drawn from the natural log, such that the natural log of the returned value is uniformly distributed. Quniform returns discrete values between the upper and lower bounds, and uniform profile returns a value uniformly between the upper and lower bounds.

The framework described in Section 3 was used to tune hyperparameters for the two load cases specified in Figure 7. The Python code developed for the framework was run on Delta supercomputer at the National Center for Supercomputing and Applications at the University of Illinois at Urbana-Champaign using an A100 single-node GPU consisting of a single AMD 64-core processor with 256 GB DDR4-3200 RAM. The optimized hyperparameters for the two load cases obtained using Tree-structured Parzen Estimator (TPE) optimization algorithms are shown in Table 4. All the cases were analyzed using 100x50 training points over the computational domain, as shown in Figure 13. Simultaneously, FEA was also undertaken using Abaqus using element type CPS4R. The results obtained from FEA were assumed to be ground truth. The deviation of displacement predicted by DEM (using tuned hyperparameters) from those obtained from FEA was calculated using the L₂-norm given by Equation 10.

$$\left|\left|L_2\right|\right|=\sqrt{\frac{\sum _{trainingpoints}\sum _{i=\mathrm{1,2}}{\left(u_{i\left(DEM\right)}-u_{i\left(FEA\right)}\right)}^2}{numberoftrainingpoints}}$$

(12)

Table 5. Optimized hyperparameters for compressive and bending loads.

Hyperparameter	Compression	Bending
Learning Rate	1.35145	1.40475
Neurons	106	98
Standard Dev (DNN)	0.01977	0.03276
Standard Dev (RFF)	0.46094	0.49815
Loss function	-4.98435	-13.34657
L2-error	0.000019	0.000037
Optimization time	8 min 13 sec	12 min 12 sec

Table 5. Optimized hyperparameters for compressive and bending loads.

Hyperparameter	Compression	Bending
Learning Rate	1.35145	1.40475
Neurons	106	98
Standard Dev (DNN)	0.01977	0.03276
Standard Dev (RFF)	0.46094	0.49815
Loss function	-4.98435	-13.34657
L2-error	0.000019	0.000037
Optimization time	8 min 13 sec	12 min 12 sec

6. Transferability of hyperparameters

Though hyperparameters can be optimized for different load cases and plate geometry, multiple DEM iterations must be undertaken, which can be time-consuming. Thus, we explored the transferability of hyperparameters to different geometries and load cases to minimize this issue. The optimized hyperparameters found during compressive loading (mentioned in Table 5) are used for the following examples.

6.1. Transferability of hyperparameters to different load cases

The optimized hyperparameters obtained for compression load case were used to predict displacements under tensile and bending loads. For the tensile test, a tensile pressure of 50 MPa was applied on the right end of the plate. Figure 14 shows predicted displacements for all three load cases. To obtain the displacements under tensile loading, 32.975 seconds were taken to train the DEM model and 0.0067 seconds to evaluate the model (to obtain displacements). Even though the parameters were optimized for compression, we notice that it can also predict deformation under uniform tensile and bending loading with L₂ errors of 1.95e-05 and 0.00192, respectively. Figure 15 shows deviations in predicted displacements for all three load cases compared to FEA.

6.2. Transferability of hyperparameters to different numerical resolution

The plate was discretized in 200 x 50 and 100 x 100 points instead of 100 x 25 points. The discretization was varied to understand the transferability of hyperparameters with variation in numerical resolution. Figure 16 shows the displacements obtained for the compressive load case in the y-direction. The figure shows that for selected hyperparameters, the displacements are not resolution dependent. To obtain the displacements for 200 x 50 and 100 x 100 points, 32.7134 seconds and 32.7906 seconds were taken to train the DEM model, and 0.0065 seconds and 0.0065 seconds were taken to evaluate the models, respectively. We notice that the resolution does not significantly affect the computational time.

6.3. Transferability of hyperparameters across geometry

A 1 mm x 1 mm plate was subjected to uniaxial compression using 25 x 25 points, as shown in Figure 17a. The training points used for DEM are shown in Figure 17b. From Figure 18, we can conclude that in DEM, the optimized hyperparameters can be used to predict displacements for geometries different from which they were optimized. 11.553 seconds were taken to train the model, while 0.0055 seconds were taken to evaluate the model. We notice that the time taken to train the model decreases as the geometry size decreases. This observation can be used in the future to reduce the time required to obtain optimal hyperparameters.

In addition, we also analyzed the transferability of DEM (with RFF mapping) to predict displacements for geometries for which it was not trained. In the following example, the DEM model was trained using a 4 mm x 1 mm plate subjected to a compressive load of -10 MPa. The trained model was then used to predict displacements produced in a 1 mm x 1 mm plate. 32.229 seconds were taken to train the model while, 0.0041 seconds were taken to evaluate the model. The results are shown in Figure 19. Even though the DEM model was trained for a plate of dimensions 4 mm x 1 mm, it can still predict displacements in the plate sized 1 mm x 1 mm subjected to similar boundary conditions.

6.4. Transferability of hyperparameters to partial loading

In this example, a load is applied locally on the upper right-hand corner of a 4 x 1 mm² plate (Figure 20). The same number of domain and boundary training points is used in Section 4. The displacement obtained from DEM, FEA, and error in prediction compared to FEA are shown in Figure 21. To obtain displacement under partial loading DEM model was trained in 43.332 seconds. 0.066 seconds were taken to obtain displacements for the model under evaluation.

Figure 21 shows that DEM can predict displacement under localized traction using the optimal hyperparameters parameters obtained from uniform compression. However, the magnitude of error increases by order of magnitude. We also notice again that the error in displacements is proportional to the magnitude of displacement in each direction.

6.5. Transferability of hyperparameters to a random geometry

In the previous section, only simplified geometries were studied. To analyze the transferability to complex geometries, passive elements were introduced by multiplying the corresponding strain energy of the element with their material density. Since passive elements do not contribute to internal energy, the material density corresponding to passive elements was zero. A random geometry of density distribution, shown in Figure 22, was introduced. A uniform compressive load of 10 MPa was applied at the right end of geometry, while the left end was fixed in all DOF. The displacements obtained from the DEM model are shown in Figure 23. The figure shows that even though the hyperparameters were optimized for a simplified rectangular plate, they can be used to predict displacements in complex geometry.

7. Effect of introducing RFF mapping

In earlier examples, RFF mapping was introduced in MLP to reduce the bias in the model’s training (as discussed in Section 3). In this section, the models are trained without RFF mapping to study its effect. As discussed in previous sections, hyperparameters are interdependent. Thus, hyperparameter tuning was undertaken for a 4x1mm² plate subjected to compressive loads. The procedure discussed in the above sections was followed. The optimized hyperparameters are shown in Table 6.

Table 6. Optimized hyperparameters for compressive load.

Hyperparameter	Compression
Learning Rate	1.00118
Neurons	98
Standard Dev (DNN)	0.05492

Table 6. Optimized hyperparameters for compressive load.

Hyperparameter	Compression
Learning Rate	1.00118
Neurons	98
Standard Dev (DNN)	0.05492

The displacements obtained from DEM without RFF mapping for three different cases: uniform compression of 4x1 mm² plate, uniform bending of 4x1 mm² plate, and uniform compression of a random geometry are shown in Figure 24, Figure 25 and Figure 26. Comparing the error in displacements under compression loading (Figure 24e, 24f) with the error when RFF mapping was present (Figure 15a, 15b), we can conclude that no significant difference exists between the two cases. However, while comparing the errors under bending and random geometry (Figure 15e, 15f with Figure 25e, 25f, and Figure 23e, 23f with Figure 26e, 26f), we conclude that displacements predicted with RFF mapping are more accurate than the displacements predicted without RFF mapping. While comparing Figure 23e, 23f with Figure 26e, 26f, we notice that the error in displacements reduces by order of magnitude. Thus, the introduction of RFF mapping improved the accuracy of DEM and enabled the transferability of hyperparameters to wider geometries and boundary conditions.

8. Conclusions

This paper proposes a two-loop architecture to obtain DEM’s best architecture and hyperparameters. The inner loop predicts displacements and potential energy using DEM based on hyperparameters chosen by an outer loop. The outer loop receives the system’s potential energy from the inner loop, which depends on the architecture of ANN defined by hyperparameters. The optimization algorithm in the outer loop then searches the non-convex optimization space to identify hyperparameters for which minimum potential energy can be achieved. Fourier feature mapping is employed in the inner loop (DEM) to improve accuracy and simplify calculations. Six hyperparameters were chosen for the study. We found that the values of hyperparameters identified using this approach and the implementation of Fourier feature mapping produce displacements similar to those obtained through FEA (with an order of accuracy of O(h^-3)).

Additionally, based on the examples presented in the paper, we conclude the following:

• The displacements obtained for tension and compression load cases are more sensitive to hyperparameters than displacements obtained for bending loads.
• With optimal hyperparameters, the order of accuracy for tension and compression load cases is lower than for bending loads.
• The optimal hyperparameters chosen through compression and tension load cases can predict displacement for different loading conditions and geometries. As a result, the optimal hyperparameters can be searched for a single loading condition and simplified geometry (depending on the desired accuracy).
• The error in predicted displacements is proportional to the magnitude of displacement.
• Generally, the activation function can be fixed to Randomized Leaky Rectified Linear Unit (rrelu) and the number of layers to 5. The two-loop architecture can then tune the hyperparameters for the given BVP.
• The accuracy of DEM and transferability of hyperparameters can be improved with RFF mapping.