Extreme sparsification of physics-augmented neural networks for
  interpretable model discovery in mechanics
Data-driven constitutive modeling with neural networks has received increased
interest in recent years due to its ability to easily incorporate physical and
mechanistic constraints and to overcome the challenging and time-consuming task
of formulating phenomenological constitutive laws that can accurately capture
the observed material response. However, even though neural network-based
constitutive laws have been shown to generalize proficiently, the generated
representations are not easily interpretable due to their high number of
trainable parameters. Sparse regression approaches exist that allow to
obtaining interpretable expressions, but the user is tasked with creating a
library of model forms which by construction limits their expressiveness to the
functional forms provided in the libraries. In this work, we propose to train
regularized physics-augmented neural network-based constitutive models
utilizing a smoothed version of $L^{0}$-regularization. This aims to maintain
the trustworthiness inherited by the physical constraints, but also enables
interpretability which has not been possible thus far on any type of machine
learning-based constitutive model where model forms were not assumed a-priory
but were actually discovered. During the training process, the network
simultaneously fits the training data and penalizes the number of active
parameters, while also ensuring constitutive constraints such as thermodynamic
consistency. We show that the method can reliably obtain interpretable and
trustworthy constitutive models for compressible and incompressible
hyperelasticity, yield functions, and hardening models for elastoplasticity,
for synthetic and experimental data.
Jan N. Fuhg, Reese E. Jones, Nikolaos Bouklas
Extreme sparsification of physics-augmented neural networks for interpretable model discovery in mechanics

###### Abstract

Data-driven constitutive modeling with neural networks has received increased interest in recent years due to its ability to easily incorporate physical and mechanistic constraints and to overcome the challenging and time-consuming task of formulating phenomenological constitutive laws that can accurately capture the observed material response. However, even though neural network-based constitutive laws have been shown to generalize proficiently, the generated representations are not easily interpretable due to their high number of trainable parameters. Sparse regression approaches exist that allow to obtaining interpretable expressions, but the user is tasked with creating a library of model forms which by construction limits their expressiveness to the functional forms provided in the libraries. In this work, we propose to train regularized physics-augmented neural network-based constitutive models utilizing a smoothed version of \(L^{0}\)-regularization. This aims to maintain the trustworthiness inherited by the physical constraints, but also enables interpretability which has not been possible thus far on any type of machine learning-based constitutive model where model forms were not assumed a-priory but were actually discovered. During the training process, the network simultaneously fits the training data and penalizes the number of active parameters, while also ensuring constitutive constraints such as thermodynamic consistency. We show that the method can reliably obtain interpretable and trustworthy constitutive models for compressible and incompressible hyperelasticity, yield functions, and hardening models for elastoplasticity, for synthetic and experimental data.

Physics-augmented machine learning Solid mechanics Data-driven constitutive models

## 1 Introduction

In continuum mechanics, material-specific constitutive models are a necessary closure relationship to describe the motion of solid bodies. In contrast to balance laws and kinematic equations, these mathematical descriptions of material behavior do not directly follow from a physical law [1]. Nevertheless, mechanistic assumptions, mathematical well-posedness, and physical understanding still constrain the formulations to adhere to objectivity, material symmetry, and thermodynamic consistency considerations, among others [2].

With the recent advances in manufacturing technologies including additive manufacturing, we are experiencing a rise of highly specific materials with advanced requirements [3]. The systematic satisfaction of these design requirements has led to a pressing need for advanced predictive capabilities, and automation directly linking experimental charac
terization and modeling. More specifically, specialized constitutive models are required that focus on constructing mathematical models capable of describing relevant physical phenomena [4] that can eventually enable computations at the structural level. In the last decades, phenomenological models have been the main driver of progress in constitutive modeling. Phenomenological models are derived and developed based on knowledge of the material response, are based on a limited number of user-chosen functional forms, and are characterized by a limited number of unknown material parameters often directly tied to specific experiments [5]. They aim to describe complex behaviors with as few parameters as possible [6] and are generally interpretable and the user has knowledge about their extrapolation behavior [7]. Extrapolation capabilities are important as traditional mechanical experiments provide stress-strain data (or rates of these quantities for rate-dependent responses) only for limited stress- or strain-states, e.g. uniaxial, biaxial, simple shear and hydrostatic, leaving the majority of the strain- or stress- space virtually unexplored. Crucially, this has been accomplished by strictly enforcing thermodynamic constraints and mechanistic assumptions in phenomenological constitutive models. However, due to the specific user-chosen model form, phenomenological models in general experience a type of model-form error [8], i.e. the model often is not descriptive enough to fully fit the data which results from incomplete knowledge about the material response [9].

In order to improve the restrictions of functional form selections of classical constitutive models and to automate their formulation, machine learning (ML)-based constitutive models have been popularized in recent years [10]. They have been used to model hyperelasticity [11, 12, 13] viscoelasticity [14, 15, 16] and plasticity [17, 18, 19]. Recently, these representations have been extended to enforce objectivity- and material symmetry-constraints [11, 20], polyconvexity [21, 22, 23, 24] and thermodynamical consistency [25, 26, 27, 28] in a physics-augmented manner. Many of these approaches are built around the flexibility of neural networks that allow them to fulfill the constitutive constraints by construction. For example, hyperelastic models as suggested by Ref. [21] are designed to be polyconvex and rely on input convex neural networks [29]; an architecture which can strictly enforce the development of a strain energy density that is polyconvex with respect to its inputs. However, even though these proposed ML-based formulations are designed to comply with mechanistic and thermodynamic constraints, they are either not easily interpretable due to the high dimensionality of parameters required in these representations (e.g. often in the order of thousands or more), or might be too restrictive if a model form (or a model-form library) is already pre-selected [30, 31, 32]. The former is also problematic in the low and limited data domain because, although the physics augmentation can act as a regularizer in some cases, the models are still overfitting as will be shown in this work. Lastly, neural network models take longer to implement in existing commercial and open-source computational infrastructure (e.g. finite element modeling platforms) that is not directly coupled to machine-learning frameworks [33]. To combat all of these points and to efficiently discover interpretable constitutive models without restrictive assumptions about specific model forms, we introduce extreme sparsification to physics-augmented neural network-based constitutive models. Our approach is based on network pruning, which is illustrated in Figure 1.

In recent years, neural network pruning - the reduction of the network size by removing parameters - has received increased interest in the machine-learning community [34, 35]. In this context, the pruning is mostly utilized to enable the networks to be deployed in real-time on mobile devices. In general, there are two ways to prune a neural network

* In a first step train a full neural network without any regularization. Then in the second step remove all trainable parameters that are below a threshold and then, in the third step, retrain the model. These greedy phases of pruning and retraining may then be repeated until the required balance between performance and network size is reached [36].
* Regularize the neural network in an eager manner directly during training without the need for any postprocessing or iterative steps [37]. This approach of course reduces the training time to obtain a pruned network but typically has a hyperparameter controlling the influence of the secondary sparsity objective.

In this work, we follow the latter since we aim to be as close as possible to the standard training procedure of physics-augmented neural networks for constitutive modeling in mechanics. In general, most techniques are built around penalizing the \(L^{p}\)-norm (\(p\geq 0\)) of the parameters using an additional term in the loss function [36, 38]. Typical norms are the \(L^{1}\)- and \(L^{2}\)-norm which are known as Lasso- and Ridge-regularization respectively [39]. Other methods prune parameters in groups to remove whole neurons or channels [40]. Recently, Ref. [41] introduced a pruning technique that aims to make neural networks more modular and interpretable by pruning the network during training by encouraging _locality_, i.e. the more neurons communicate the closer they should be in Euclidean space. This is achieved by a local \(L^{1}\)-regularization measure and by placing neurons with high communication closer together by employing a swapping algorithm. In this work, we rely on Ref. [42] which prune the network through a smoothed version of the expected \(L^{0}\)-regularization. The benefits of this approach include that it enforces sparsity without placing a penalty on the magnitude of the weights and it allows for parameters to be exactly zero. Hence, no thresholding is necessary. We remark that our work is related to neural symbolic regression [43, 44], we, however, aim to use standard neural network models that have been used in the data-driven constitutive modeling community and combine these formulations with
sparsification techniques instead of selecting a model from a library of formulations. We choose this approach to enhance the expressiveness of the models and avoid selecting a model-informed basis for our representation.

The paper is structured as follows. Section 2 introduces physics-augmented neural networks and explains their sparsification through smoothed \(L^{0}\) regularization. We test the framework's ability to produce interpretable constitutive models for hyperelasticity and elastoplasticity in Section 3. The paper is concluded in Section 4.

## 2 Sparsifying physics-augmented neural networks for constitutive modeling

Due to their remarkable flexibility, neural networks have emerged as the most used machine-learning approach for data-driven constitutive modeling. In particular, compared to other regression techniques, such as Gaussian process regression or Support Vector Regression that have been applied for constitutive modeling [45, 46, 6], neural networks have received more attention with regards to finding a proficient balance between enforcing constitutive-constraints and retaining model expressiveness [47]. In the following, we will show two examples of incorporating constraints into neural networks and discuss a potential way of sparsifying them. We remark that in this context we mean the sparsification of all trainable parameters of a neural network, instead of, e.g., the sparse identification of terms from a user-defined functional library as in SINDy [48, 49] or the EUCLID framework [7, 50].

### Physics-augmented neural network formulations

Two classes of physics augmentation have been found to be crucial for constitutive modeling

* Input convex/concave functions. Convex functions are needed to model polyconvex strain energy density functions for hyperelasticity [21, 24] and guarantee dissipation requirements when utilizing potentials to model yield functions [51, 19, 52] or hardening behavior [53, 27].

Figure 1: Extreme sparsification of a physics-augmented neural network constitutive model to improve their interpretability. Note the final input convex network utilizes the pass-through of \(I_{1}\) and \(I_{2}\) terms.
* Positive, monotonically increasing functions or their counterparts negative, monotonically decreasing functions. These functions are needed to model the derivatives of convex functions [26] or where mechanistic assumptions require monotonic effects such as isotropic hardening [27].

Even though these examples of constraints might sound limited, they form the basis of the framework for constitutive modeling for several classes of materials, whose responses are described via hyperelasticity, viscoelasticity, elastoplasticity, viscoplasticity and damage mechanics [54]. We will briefly highlight how both of these constraints can be intrinsically incorporated into feed-forward neural networks.

#### 2.1.1 Input convex neural network

Following Ref. [29], consider an output \(\hat{\mathbf{y}}\in\mathbb{R}^{n^{L}}\) connected to an input \(\mathbf{x}_{0}\in\mathbb{R}^{n^{0}}\) by the neural network \(\mathcal{N}\) given by

\[\mathbf{x}_{1}=\sigma_{1}\left(\mathbf{x}_{0}\mathbf{W}_{1}^{T}+\mathbf{b}_{1}\right) \in\mathbb{R}^{n^{1}} \tag{1}\] \[\mathbf{x}_{l}=\sigma_{l}\left(\mathbf{x}_{l-1}\mathbf{W}_{1}^{T}+\mathbf{x}_{0} \mathbf{\mathcal{W}}_{1}^{T}+\mathbf{b}_{l}\right) \in\mathbb{R}^{n^{l}},\qquad l=2,\ldots,L-1\] \[\hat{\mathbf{y}}=\mathbf{x}_{L-1}\mathbf{W}_{L}^{T}+\mathbf{x}_{0}\mathbf{\mathcal{W} }_{L}^{T}+\mathbf{b}_{L}, \in\mathbb{R}^{n^{L}}\]

with the weights \(\mathbf{W}\) and \(\mathbf{\mathcal{W}}\), the biases \(\mathbf{b}\) and the activation functions \(\sigma\). The weights and biases form the set of trainable parameters \(\mathbf{\theta}=\{\{\mathbf{W}_{i}\}_{i=1}^{L},\{\mathbf{\mathcal{W}}_{i}\}_{i=2}^{L},\{ \mathbf{b}_{i}\}_{i=1}^{L}\}\). The output is then convex with regards to the input if the weights \(\{\mathbf{W}_{i}\}_{i=2}^{L}\) are non-negative and the activation functions \(\{\sigma_{i}\}_{i=1}^{L}\) are convex and non-decreasing. For a proof see Ref. [29].

#### 2.1.2 Positive, monotonically increasing neural network

Let a positive, monotonically increasing neural network \(\mathcal{N}\) be defined by

\[\mathbf{x}_{0} \in\mathbb{R}^{n^{0}} \tag{2}\] \[\mathbf{x}_{1}=\sigma_{1}\left(\mathbf{x}_{0}\mathbf{W}_{1}^{T}+\mathbf{b}_{1}\right) \in\mathbb{R}^{n^{1}}\] \[\mathbf{x}_{l}=\sigma_{l}\left(\mathbf{x}_{l-1}\mathbf{W}_{1}^{T}+\mathbf{b}_{l}\right) \in\mathbb{R}^{n^{l}},\qquad l=2,\ldots,L-1\] \[\hat{\mathbf{y}}=\mathbf{x}_{L-1}\mathbf{W}_{L}^{T}+\mathbf{b}_{L}, \in\mathbb{R}^{n^{L}}\]

where \(\mathbf{W}\), \(\mathbf{b}\) and \(\sigma\) denote the weights, biases and activation functions respectively. The set of trainable parameters reads \(\theta=\{\{\mathbf{W}_{i=1}^{L}\},\{\mathbf{b}_{i}\}_{i=1}^{L}\}\). Then each output value of \(\hat{\mathbf{y}}\) is positive, and monotonically increasing with regard to all input values of \(\mathbf{x}_{0}\) when the trainable parameters are nonnegative and the activation functions are positive and non-decreasing. For a proof, we refer to Ref. [27].

### Extreme sparsification with smoothed \(L^{0}\) regularization

Given a data set of input-output pairs \(\{\mathbf{x}^{i},\mathbf{y}^{i}\}_{i=1}^{N}\) the trainable parameters of eqs. (1) and (2) are classically found by minimizing a loss function \(\mathcal{R}(\mathbf{\theta})\)

\[\mathbf{\theta}^{\star}=\operatorname*{arg\,min}_{\mathbf{\theta}}\mathcal{R}(\mathbf{ \theta})=\operatorname*{arg\,min}_{\mathbf{\theta}}\frac{1}{N}\sum_{i=1}^{N}\left[ \mathcal{L}(\mathcal{N}(\mathbf{x}^{i};\mathbf{\theta}),\mathbf{y}^{i})\right]\]

with some loss function \(\mathcal{L}(\bullet)\). In the following, we discuss how to sparsify the parameters by adding regularization. The approach follows the general idea of a gating system where each trainable parameter is multiplied by a gate value \(z\in[0,1]\) of either zero or one depending on whether the parameter should be active or not. Zero reflects an inactive parameter while one defines an active parameter. The number of active gates could then be penalized in the loss function. However, due to the binary nature of the gates, the loss function would not be differentiable. Hence, following Ref. [42] we consider a reparametrization of the trainable parameters using a smoothed "gating" system, i.e. let

\[\mathbf{\theta}=\overline{\mathbf{\theta}}\odot\mathbf{z},\quad\text{with}\quad\mathbf{z}= \min(\mathbf{1},\max(\mathbf{0},\overline{\mathbf{s}})) \tag{3}\]

where \(\odot\) denotes the Hadamard product and where

\[\overline{\mathbf{s}}=\mathbf{s}(\zeta-\gamma)+\gamma,\quad\mathbf{s}=\text{Sigmoid}(( \log\mathbf{u}-\log(1-\mathbf{u})+\log\mathbf{\alpha})/\beta),\quad\mathbf{u}\sim U(\mathbf{0}, \mathbf{1}). \tag{4}\]

Here, \(\gamma\), \(\beta\), \(\zeta\) and \(\log\mathbf{\alpha}\) are user-chosen parameters that define the smoothing of the gate \(\mathbf{z}\). Following the suggestions of the authors [42], we choose \(\gamma=-0.1\), \(\zeta=1.1\), \(\beta=2/3\) and \(\log\mathbf{\alpha}\) is obtained by sampling from a normal distribution with zero mean and a standard deviation of \(0.01\).
Since the gate values are dependent on a random variable, we can define a Monte Carlo approximated loss function as

\[\mathcal{R}(\overline{\mathbf{\theta}})=\frac{1}{L}\sum_{l=1}^{L}\left(\frac{1}{N} \left(\sum_{i=1}^{N}\mathcal{L}(\mathcal{N}(\mathbf{x}_{i},\overline{\mathbf{\theta}} \odot\mathbf{z}^{l}),\mathbf{y}_{i})\right)\right]+\lambda\sum_{j=1}^{|\theta|}\text{ Sigmoid}(\log\alpha_{j}-\beta\log\frac{-\gamma}{\zeta}) \tag{5}\]

with \(L\) being the number of samples of the Monte Carlo approximation and where \(\lambda\) is a weighting factor for the regularization. At test time we can then set the values of the trainable parameters \(\mathbf{\theta}^{\star}=\overline{\mathbf{\theta}}^{\star}\odot\hat{\mathbf{z}}\) where the gates can be obtained as

\[\hat{\mathbf{z}}=\min(\mathbf{1},\max(\mathbf{0},\text{Sigmoid}(\log\mathbf{\alpha})(\zeta- \gamma)+\gamma)). \tag{6}\]

In the following section, we show how this idea can be applied to constitutive modeling in computational mechanics.

## 3 Results

We consider three different application areas where fitting physics-augmented neural network constitutive laws have been applied in the literature:

* Incompressible and compressible hyperelasticity
* Yield functions for elastoplasticity
* Isotropic hardening for elastoplasticity.

For each case, we look at multiple different experimental and synthetic datasets to test the capabilities of the sparsification approach. We have implemented all the models in Pytorch [55] and use the Adam optimizer [56] with a learning rate of \(10^{-3}\) to solve the optimization problem. We also note that for elastoplasticity, we follow the paradigm suggested in [27] for modular learning approaches.

### Constitutive modeling of hyperelastic potential

Hyperelasticity enables modeling finite strain, nonlinear elastic material behavior by assuming that a strain-dependent potential function \(\Psi(\mathbf{F})\) exists from which the stress can be derived [2]. Here \(\mathbf{F}=\frac{\partial\phi}{\partial\mathbf{X}}\) is the deformation gradient that depends on \(\phi(\mathbf{X},t)\) which describes the motion of a body from a reference position \(\mathbf{X}\) to its current position at time \(t\). One important strain measure for finite strain modeling is the right Cauchy-Green tensor \(\mathbf{C}=\mathbf{F}^{T}\mathbf{F}\). In the following, we restrict ourselves to isotropic material behavior. The hyperelastic potential is generally subject to constitutive constraints that are derived from thermodynamic considerations and mechanistic assumptions [21]. These include

* Objectivity and material symmetry. Let \(\mathbf{R}\) be an orthogonal matrix, the strain energy function then needs to adhere to the assumed isotropic nature of the material behavior, i.e. \[\Psi(\mathbf{F}\mathbf{R}^{T})=\Psi(\mathbf{F}).\] (7) Furthermore, the energy value needs to be independent of the choice of an observer, which is known as objectivity, and which is defined by \[\Psi(\mathbf{R}\mathbf{F})=\Psi(\mathbf{F}).\] (8) By formulating the strain energy function in terms of the invariants of the right Cauchy-Green tensor \(\psi(I_{1},I_{2},J)\) where \[I_{1}=\text{tr}(\mathbf{C}),\quad I_{2}=\text{tr}(\text{cof}\mathbf{C}),\quad J=\sqrt {\det(\mathbf{C})},\] (9) instead of the deformation gradient, the objectivity and the material symmetry conditions are fulfilled.
* Normalization condition. It is physically sensible to define the form of the strain energy function such that both its value and the derived stress value are zero at the undeformed configuration \(\mathbf{C}=\mathbf{I}\)[2].
* Polyconvexity. Following Ball [57], polyconvex hyperelasticity guarantees the existence of minimizers of the underlying potential in finite elasticity. Polyconvexity in this context requires that \(\Psi(\mathbf{F},\text{cof}\mathbf{F},\det\mathbf{F})\) is convex in all its arguments [58]. This is equivalent to ensuring convexity of the strain energy function with regard to a set of polyconvex invariants.

In the following, we discuss neural network formulations that enforce all these constitutive conditions and use regularization to obtain interpretable model forms for different data sets. We differentiate between problems in compressible and incompressible hyperelasticity.
#### 3.1.1 Compressible hyperelastic materials from numerical data

In order to generate a compressible, polyconvex, objective, and isotropic material description the strain energy function can be assumed to be a convex function \(\psi(I_{1},I_{2},J)\). We then build an input convex neural network, c.f. eq. (1) that takes these invariants as input and outputs a scalar \(\hat{\Psi}^{NN}(I_{1},I_{2},J)\). In order to also fulfill the normalization condition, the approach suggested by Ref. [24] can be used, we thus obtain the strain energy prediction as

\[\hat{\Psi}(I_{1},I_{2},J)=\hat{\Psi}^{NN}(I_{1},I_{2},J)-\hat{\Psi}^{NN}(3,3,1)- \psi^{S}(J) \tag{10}\]

where \(\psi^{S}\) defined as

\[\psi^{S}=n(J-1)=\underbrace{\left.\left(2\frac{\partial\hat{\Psi}^{NN}}{ \partial I_{1}}+4\frac{\partial\hat{\Psi}^{NN}}{\partial I_{2}}+\frac{ \partial\hat{\Psi}^{NN}}{\partial J}\right)\right|_{\mathbf{C}=\mathbf{I}}}_{=n}(J-1) \tag{11}\]

enforces the normalization of the stress. For more information see Appendix A.1. A prediction of the second Piola-Kirchhoff stress can then be obtained with

\[\hat{\mathbf{S}}=2\frac{\partial\hat{\Psi}}{\partial\mathbf{C}}=2\left(\sum_{i}\frac{ \partial\hat{\Psi}}{\partial I_{i}}\frac{\partial I_{i}}{\partial\mathbf{C}} \right)=2\left(\frac{\partial\hat{\Psi}}{\partial I_{1}}+I_{1}\frac{\partial \hat{\Psi}}{\partial I_{2}}\right)\mathbf{I}-2\frac{\partial\hat{\Psi}}{\partial I _{2}}\mathbf{C}+J\frac{\partial\hat{\Psi}}{\partial J}\mathbf{C}^{-1}. \tag{12}\]

Using this formulation we can train a sparsified, compressible, physics-augmented neural network model on stress-strain data. For this application, consider a region in the deformation gradient space of the form

\[F_{ij}\in[F_{ij}^{L},F_{ij}^{U}]=\begin{cases}[1-\delta,1+\delta]&\text{if}\,i =j\\ [-\delta,+\delta]&\text{else}.\end{cases} \tag{13}\]

Following Ref. [59] we define a training region with \(\delta=0.2\) and a test region with \(\delta=0.3\). We then utilize the space-filling sampling algorithm proposed in Ref. [20] to sample \(50\) training and \(10,000\) testing data points in their respective regions. In terms of the invariants, the 50 data points are summarized in Table A1. We furthermore divide the training dataset in an \(80/20\) split into training and validation points.

ExamplesFirstly, we study the performance of the proposed approach and how it is influenced by the network size and the weighting factor \(\lambda\). Consider a Gent-Gent hyperelastic model [60] of the form

\[\Psi_{gent}=-\frac{\theta_{1}}{2}J_{m}\log\left(1-\frac{I_{1}-3}{J_{m}} \right)-\theta_{2}\log\left(\frac{I_{2}}{J}\right)+\theta_{3}\left(\frac{1}{2 }(J^{2}-1)-\log J\right) \tag{14}\]

where we employ \(\theta_{1}=2.4195\), \(J_{m}=77.931\), \(\theta_{2}=-0.75\) and \(\theta_{3}=1.20975\), see Ref. [59]. We consider three different architectures: 1 hidden layer with 30 neurons, 2 hidden layers with 30 neurons, and 3 hidden layers with 30 neurons), 5 different regularization parameter values, and repeat the training process \(10\) times with different random seeds. Figure 2 depicts the average (in solid lines) and each individual data point of the final training loss and the number of parameters with an absolute value greater than \(0\) at the end of \(100,000\) training epochs for each regularization parameter value and architectural setting. We additionally highlight the average training loss (red dotted line) and the number of active parameters (in the legend) when no regularization was employed. We can see that the final training loss value as well as the number of active training parameters is relatively similar in terms of absolute values for all architectures. As was expected, higher regularization parameters penalize the number of active parameters more, and hence we see can see that the number of active parameters increases with lower regularization parameters. Since the loss is then also dominated by the regularization, the training loss decreases with a decrease in \(\lambda\). Due to the higher number of tunable parameters in the unregularized models (e.g. \(1112\) parameters whose absolute value is greater than zero for the architecture with 1 hidden layer and 30 neurons), these networks achieve a comparatively low training loss. Even though, the training loss for the sparsified networks with \(\lambda=10^{-5}\) might reach a similar value with around \(15\) parameters.

For the same trained networks as above, Figure 3 summarizes the generalization performance of the networks based on the test loss over the \(10,000\) test points. Here, again, the red dotted lines indicate the average test loss when no regularization was employed. Interestingly, the average test loss of the networks with a regularization value of \(\lambda=10^{-5}\) is better than that of the full models where \(\lambda=0\). This suggests that due to the additional regularization, these models are not only more interpretable but they also generalize more proficiently. The functional form of the network (1 hidden layer and \(\lambda=10^{-4}\)) with the median training loss is given by

\[\hat{\Psi}=0.398J+3.095\log\left(\left(1+e^{-1.356I_{2}}\right)^{1.314}\left(e^ {0.755I_{1}}+1\right)^{0.515}\left(e^{0.135I_{1}-0.319I_{2}-0.329J}+1\right)^{1. 874}+1\right)-6.686. \tag{15}\]
The representation consists of 3 terms that can be written out in one line of text, making it easy to implement and allowing for easier interpretation than the equivalent \(1112\) parameters of the full model. For example, it is easy to see that the obtained hyperelastic law does not fulfill the coercivity condition [21], i.e., as \(J\to 0^{+}\), we do not get \(\hat{\Psi}\rightarrow\infty\). This observation would not have been possible with over \(1000\) parameters.

We remark that we could have included a functional term in the strain energy function, such that the coercivity would be fulfilled by design, c.f. Ref. [21].

To highlight how accurate the obtained formulation nevertheless is we can look at the stress-strain behavior inside and outside of the training regime \(0.6\leq F_{11}\leq 1.4\), see Figure 4, where the green lines indicate the training domain. We can see that the network finds a sparse representation that appears to interpolate well and is also proficient in extrapolation.

Next, we look at the effect of the physics augmentation compared to a standard neural network model. Figure 5 depicts the influence of the \(\mathcal{L}^{0}\) regularization parameter and the network architecture on the averaged (over 10 runs) test loss and number of active parameters of the Gent-Gent model. Two things are noteworthy: (i) the test loss is generally similar between the two network formulations, and (ii) the average number of active parameters of the standard neural networks seems to be significantly more sensitive to the regularization parameter than its physics-augmented counterparts. Apart from the guaranteed adherence to physical principles and mechanistic assumptions, this seems to be a major positive side effect of the physics augmentation. The robustness of the physics-augmented approach to maintain interpretability irrespective of the choice of the regularization parameter is a key feature of the proposed framework here, and it will positively impact the overall trustworthiness of the approach, as will be discussed further in the manuscript.

Lastly, to prove the flexibility of the proposed approach, we also consider two more hyperelastic laws, in particular, a compressible Mooney-Rivlin model [61] and a polynomial model [7]. We again ran 10 random realizations of the training process and report the result with the median final training loss after \(100,000\) epochs. The true functional

Figure 3: Gent-Gent law test loss for different architectures. The red dotted line indicates the test loss of the network when no \(\mathcal{L}^{0}\)-sparsification is employed.

Figure 2: Gent-Gent law training loss for different architectures. The red dotted line indicates the training loss of the network when no \(\mathcal{L}^{0}\)-sparsification is employed.
forms as well as their sparsified representations of the neural network trained on 40 training data points are shown in Table 1. We can see that these fitted functional forms are similarly interpretable to the one discussed for the Gent-Gent model. For example, neither of the other two obtained representations fulfill the coercivity condition either.

Furthermore, similarly, to the Gent-Gent representation, the formulations obtained by these networks also appear accurate inside and even far outside the training regime, see Figure 6, where the training domain is highlighted by the green dotted lines.

### Incompressible hyperelastic materials from experimental data

If the material is incompressible, it can accommodate no volumetric deformations, captured by the constraint \(J=1\). We can then assume a predictor of the strain energy function given by

\[\hat{\Psi}(I_{1},I_{2})=\hat{\Psi}^{NN}(I_{1},I_{2})-(p+n)(J-1)-\hat{\Psi}^{NN }(3,3) \tag{16}\]

where the second term on the right-hand side of this expression has no contribution on the strain energy, \(p\) is the hydrostatic pressure (remaining to be determined by the boundary value problem at hand) and

\[n=2\left.\left(\frac{\partial\hat{\Psi}^{NN}}{\partial I_{1}}+2\frac{\partial \hat{\Psi}^{NN}}{\partial I_{2}}\right)\right|_{\mathcal{C}=\mathbf{I}} \tag{17}\]

Figure 4: Stress strain curves of fitted Gent-Gent hyperelastic laws. Green dotted lines represent the limits of the training regime.

Figure 5: Effect of the network architecture and the \(\mathcal{L}^{0}\)-regularization parameter on the Gent-Gent law test loss and the number of active parameters for physics-augmented neural networks and standard neural networks. Averaged over 10 runs.
**Gent-Gent:**

\[\begin{array}{c}\Psi_{gent}=-\frac{\theta_{1}}{2}J_{m}\log\left(1-\frac{I_{1}-3 }{J_{m}}\right)-\theta_{2}\log\left(\frac{I_{2}}{J}\right)+\theta_{3}\left( \frac{1}{2}(J^{2}-1)-\log J\right),\\ \theta_{1}=2.4195,J_{m}=77.931,\theta_{2}=-0.75,\theta_{3}=1.20975\end{array}\]

**Gent-Gent fit:**

\[\hat{\Psi}=0.398J+3.095\log\left(\left(1+e^{-1.356I_{2}}\right)^{1.314}\left(e ^{0.755I_{1}}+1\right)^{0.515}\left(e^{0.135I_{1}-0.319I_{2}-0.329J}+1\right)^ {1.874}+1\right)-6.686\]

**Mooney-Rivlin:**

\[\begin{array}{c}\Psi_{MR}=\theta_{1}\left(\frac{I_{1}}{J^{2/3}}-3\right)+ \theta_{2}\left(\frac{I_{2}}{J^{3/3}}-3\right)+\theta_{3}(J-1)^{2}\;,\\ \theta_{1}=9.2\cdot 10^{-4},\theta_{2}=2.37\cdot 10^{-3},\theta_{3}=10.0010\end{array}\]

**Mooney-Rivlin fit:**

\[\hat{\Psi}_{MR}=90.338J+17.263\log\left(\left(1+e^{-0.368J}\right)^{34.796}+1 \right)+0.276\log\left(\left(e^{-1.671I_{1}+1.305I_{2}}+1\right)^{1.298}+1 \right)-406.511\]

**Polynomial:**

\[\begin{array}{c}\Psi_{poly}=\theta_{1}\left(I_{1}-3\right)^{2}+\theta_{2} \left(I_{1}-3\right)^{4}+\theta_{3}\left(I_{2}-3\right)^{2}+\theta_{4}\left(I _{2}-3\right)^{4}+\theta_{5}\left(I_{3}-1\right)^{2},\\ \theta_{1}=0.1,\theta_{2}=0.15,\theta_{3}=2\cdot 10^{-4},\theta_{4}=1\cdot 10^{-4}, \theta_{5}=0.125\end{array}\]

**Polynomial fit:**

\[\begin{array}{c}\hat{\Psi}_{poly}=-0.116J+9.046\log\left(\left(1+e^{-0.179I _{2}}\right)^{3.59}\left(e^{0.193I_{2}}+1\right)^{1.927}+1\right)\\ \hskip 142.26378pt+0.597\log\left(\left(e^{-2.931I_{1}+2.272I_{2}}+1\right)^{ 0.696}+1\right)-33.35\end{array}\]
is determined by enforcing the normalization constraint. Furthermore \(\hat{\Psi}^{NN}(\bullet)\) is the output of an input convex neural network. Formulation (16) allows us to obtain a model that fulfills all constitutive constraints described above. The derivation can be found in Appendix A.2.

We can then obtain a formulation for the first Piola-Kirchoff stress with

\[\mathbf{P}=2\left(\left[\frac{\partial\hat{\Psi}}{\partial I_{1}}+I_{1}\frac{ \partial\hat{\Psi}}{\partial I_{2}}\right]\mathbf{F}-\frac{\partial\hat{\Psi}}{ \partial I_{2}}\mathbf{F}\mathbf{C}\right)-(p+n)J\mathbf{F}^{-T} \tag{18}\]

which in terms of the principal strains \(\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\) reads

\[P_{i}\lambda_{i}=2\left(\lambda_{i}^{2}\frac{\partial\hat{\Psi}}{\partial I_{ 1}}+\frac{\partial\hat{\Psi}}{\partial I_{2}}\left[(\lambda_{1}^{2}+\lambda_{ 2}^{2}+\lambda_{3}^{2})\lambda_{i}^{2}-\lambda_{i}^{4}\right]\right)-p-n. \tag{19}\]

Analytical forms for the hydrostatic pressure can be found for some simple boundary value problems.

In the following, we fit curves to data from uniaxial tension (UT), equibiaxial tension (ET), pure shear stress (PS), simple shear deformation (SS), and simple torsion (ST). The relevant deformations and quantities of interest that were used to fit the data are summarised in Table 2 and described in detail in Appendix A.2.

### Examples

In the following, we analyze and study the performance of the sparsified physics-augmented ML model on different sets of experimental data. We have repeated the training process 10 times with different random seeds and in the following only report the results corresponding to the median training loss after (depending on the dataset) \(50,000\) or \(80,000\) epochs. We utilize an architecture of 1 hidden layer with 30 neurons and \(\lambda=10^{-3}\) for all the following results.

\begin{table}
\begin{tabular}{|c||c|} \hline Type & Deformation & Relevant experimental output \\ \hline \hline Uniaxial tension (UT) & \(\mathbf{F}=\begin{bmatrix}\lambda_{1}&0&0\\ 0&\frac{1}{\sqrt{\lambda_{1}}}&0\\ 0&0&\frac{1}{\sqrt{\lambda_{1}}}\end{bmatrix}\) & \(P_{1}=2\left(\frac{\partial\hat{\Psi}}{\partial I_{1}}+\frac{1}{\lambda_{1}} \frac{\partial\hat{\Psi}}{\partial I_{2}}\right)\left[\lambda_{1}-\frac{1}{ \lambda_{1}^{4}}\right]\) \\ \hline Equibiaxial tension (ET) & \(\mathbf{F}=\begin{bmatrix}\lambda_{1}&0&0\\ 0&\lambda_{1}&0\\ 0&0&\frac{1}{\lambda_{1}^{2}}\end{bmatrix}\) & \(P_{1}=P_{2}=2\left(\frac{\partial\hat{\Psi}}{\partial I_{1}}+\lambda_{1}^{2} \frac{\partial\hat{\Psi}}{\partial I_{2}}\right)\left[\lambda_{1}-\frac{1}{ \lambda_{1}^{4}}\right]\) \\ \hline Pure shear stress (PS) & \(\mathbf{F}=\begin{bmatrix}\lambda_{1}&0&0\\ 0&1&0\\ 0&0&\frac{1}{\lambda_{1}}\end{bmatrix}\) & \(P_{1}=2\left(\frac{\partial\hat{\Psi}}{\partial I_{1}}+\frac{\partial\hat{\Psi} }{\partial I_{2}}\right)\left[\lambda_{1}-\frac{1}{\lambda_{1}^{4}}\right]\) \\ \hline Simple shear deformation (SS) & \(\mathbf{F}=\begin{bmatrix}1&\gamma&0\\ 0&1&0\\ 0&0&1\end{bmatrix}\) & \(P_{12}=2\gamma\left(\frac{\partial\hat{\Psi}}{\partial I_{1}}+\frac{\partial \hat{\Psi}}{\partial I_{2}}\right)\) \\ \hline Simple torsion & \(\mathbf{F}=\begin{bmatrix}1&0&0\\ 0&1&\rho\phi\\ 0&0&1\end{bmatrix}\) & \(\tau=\int_{0}^{1}4\pi\rho^{3}\phi\left(\frac{\partial\hat{\Psi}}{\partial I_{ 1}}+\frac{\partial\hat{\Psi}}{\partial I_{2}}\right)d\rho\) \\ \hline \end{tabular}
\end{table}
Table 2: Studied deformation modes with relevant output
#### 3.3.1 Model for rubber

One of the most used datasets for incompressible hyperelasticity goes back to Treloar [62] who reported the uniaxial, equibiaxial and pure shear responses of vulcanized rubber for \(20^{\circ}C\) and \(50^{\circ}C\). Following Ref. [63] this data is summarized in Table A2. In the following, we use the UT and ET data as training data and test the performance on the PS dataset.

\(20^{\circ}C\) datasetFor the \(20^{\circ}C\) dataset, the evolution of the median training loss and the median number of active parameters is shown in Figure (a)a. It can be seen that both values reduce simultaneously and that the number of active parameters reduces from an initial value greater than \(1000\) to around \(10\) over the training process. The training data, the predicted stress responses, and the \(R^{2}\)-scores are depicted in Figure (b)b. The predicted responses on the training data are proficiently accurate. Further, even on the unseen PS data the model reaches an accuracy \(R^{2}>0.99\). The obtained sparsified representation for the compressible hyperelastic law reads

\[\hat{\Psi}=0.15I_{1}+6.3\cdot 10^{-3}I_{2}-0.33J-0.28-p\left(J-1.0\right)+794. 94\log\left(10^{-4}\left(e^{0.09I_{1}-0.01I_{2}}+1\right)^{0.73}+1\right) \tag{20}\]

which fulfills all the discussed constraints such as convexity and the normalization conditions.

\(50^{\circ}C\) datasetFor the \(50^{\circ}C\) dataset, the evolution of the training loss and the number of active parameters is plotted in Figure (a)a. At the end of the training process, the discovered strain formulation is given by

\[\hat{\Psi} =0.2025I_{1}+0.0114I_{2}-0.4512J-p\left(J-1.0\right)+5.1533\log \left(0.0009\left(0.2155e^{0.1852I_{1}-0.0028I_{2}}+1\right)^{0.7778}+1\right) \tag{21}\] \[+12.5236\log\left(0.0009e^{0.0052I_{2}}+1\right)-0.2076.\]

Using this formulation, Figure (b)b shows the stress responses, the training data points, and the \(R^{2}\)-scores. Here again, the predicted responses are accurate and the model shows it can generalize well. For the same dataset Figure 9 depicts the effect of the regularization. We can see that if we train a physics-augmented neural network with 1 hidden layer with \(30\) neurons without reducing the number of trainable parameters with the \(\mathcal{L}^{0}\)-regularizer, the number of trainable parameters of the run with the median final training loss (over 10 runs) is still in the hundreds at the end of the training process. Furthermore, while the training data is fitted accurately in a comparable manner to the obtained model with \(\mathcal{L}^{0}\)-regularization, the model generalizes significantly worse on the unseen PS data, c.f. Figure (b)b. We can hence see that in the limited data domain, the regularization not only increases the interpretability of the model but improves its generalization performance as well.

Figure 7: Rubber response \(20^{\circ}C\). (a) Training loss and number of active parameters over the training process, (b) Data fit and \(R^{2}\) errors.
#### 3.3.2 Model for human brain tissue

Next, we employ the \(\mathcal{L}^{0}\)-regularized neural network approach to find incompressible hyperelastic formulations for human brain tissues. In particular, we look at four datasets of mechanical tests on the cortex, the corona radiata and two different sections of the midbrain.
we argue that its complexity is (roughly) in line with other reported forms used to model human brain tissues, e.g. compare Ref. [66] which use a 3-term Ogden model [67]. Figure 9(b) depicts the data points, the predicted responses outside, and the \(R^{2}\)-error values. We can see that not only the training data is fitted proficiently well, and the model is also highly accurate on the shear stress test data.

The loss behavior and the number of active parameters of the network fitted on the Corona Radiata are shown in Figure 10(a). We can see a similar performance to the Cortex example. The final functional form reads

\[\begin{split}\hat{\Psi}_{\text{CR Radiata}}&=-43.126J-p \left(J-1.0\right)-56.327\\ &+4962.449\log\Big{(}0.014\left(1+0.0001e^{-1.294I_{1}}\right)^{ 66688.773}\left(0.0001e^{0.841I_{2}}+1\right)^{189.231}+1\Big{)}.\end{split} \tag{23}\]

The fitted responses and the \(R^{2}\)-score are shown in Figure 10(b) where again proficient generalization of the obtained model can be noted.

MidbrainNext, we fit a model to UT and ST data for two sections of the midbrain. The datsets were adopted from Ref. [31] and are listed in Table A4. For the first section, Figure 11(a) shows the training loss corresponding to the

Figure 11: Corona radiata dataset. (a) Training loss and number of active parameters over the training process, (b) Data fit and \(R^{2}\) errors.

Figure 10: Cortex dataset. (a) Training loss and number of active parameters over the training process, (b) Data fit and \(R^{2}\) errors.
training run with the median final training loss and the respective evolution of the number of active parameters. The obtained model form for this section of the midbrain reads

\[\begin{split}\Psi_{\text{MB1}}&=-47.0604J-p\left(J-1. 0\right)+28.8387\\ &+9.9229\log\Big{(}0.0048\big{(}1+0.0001e^{-0.4491I_{1}}\big{)}^{ 194016.0037}\left(0.0001e^{1.299I_{2}}+1\right)^{399.3928}+1\Big{)}.\end{split} \tag{24}\]

We remark that due to the lower number of parameters the loss appears to have single optimization steps where the training loss sharply increases but then immediately reverts. This could potentially be avoided by a lower learning rate, which is however out of the scope of the current paper, and appears to have no effect on the final accuracy of the obtained model. The latter is highlighted in Figures 11(b) and 11(c) that show the data and the model response for UT and ST, respectively. We furthermore highlight the \(R^{2}\) score which for both cases is \(>0.97\) and depict the fit of Flaschel et. al. [31] which was obtained using a sparse regression framework called EUCLID [7]. Ref. [31] obtain a strain energy function of the form

\[\Psi_{\text{Fl}}=\frac{2\cdot 0.01}{(-78.58)^{2}}\left(\lambda_{1}^{-78.58}+ \lambda_{2}^{-78.58}+\lambda_{3}^{-78.58}-3\right)+\frac{2\cdot 90.36}{(-28.71)^{2} }\left(\lambda_{1}^{-28.71}+\lambda_{2}^{-28.71}+\lambda_{3}^{-28.71}-3\right) \tag{25}\]

which could potentially be perceived as slightly more interpretable than the strain energy function of eq. (24), which however has a significantly worse accuracy with \(R^{2}\)-scores below \(0.86\). The loss and parameter evolution of the

Figure 12: Midbrain Section 1 response. (a) Training loss and number of active parameters over the training process, (b) Uniaxial tension data, the fit of the proposed approach, and fit of Ref. [31], (c) torsion data, the fit of the sparsified neural network approach, compared to the fit of Ref. [31].
network trained on the data of the second midbrain section is depicted in Figure 12(a). The final model form reads

\[\begin{split}\Psi_{\text{MB2}}&=-57.0237J-p\left(J-1.0 \right)+29.8625\\ &+87.1968\log\left(0.0015\left(1+0.0003e^{-0.2335I_{1}}\right)^{ 24578.9446}\left(0.0003e^{0.7607I_{2}}+1\right)^{437.4412}+1\right)\end{split} \tag{26}\]

which is of similar complexity to the model obtained for the first midbrain section. Figures 12(b) and 12(c) plot the data points, our fit and the fit of Ref. [31]. Again the proposed approach shows a proficient accuracy which is better than the one provided by the EUCLID method.

### Constitutive modeling of elastoplastic material responses

Elastoplasticity is a framework to model history-dependent nonlinear materials. In the small strain regime, we assume that we can split the strain into elastic and plastic strains \(\mathbf{\epsilon}=\mathbf{\epsilon}^{e}+\mathbf{\epsilon}^{p}\). We can then postulate the existence of a free energy function \(\Psi(\mathbf{\epsilon}^{e},r)=\Psi^{e}(\mathbf{\epsilon}^{e})+\Psi^{r}(r)\) that is decomposed into its elastic and part \(\Psi^{e}\) and a history-dependent component that models isotropic hardening \(\Psi^{r}\) which is dependent on the internal variable \(r\). To be consistent with thermodynamics, the intrinsic dissipation has to be non-negative [68], i.e.

\[\mathcal{D}_{int}=\mathbf{\sigma}:\dot{\mathbf{\epsilon}}-\dot{\Psi}\geq 0= \left(\mathbf{\sigma}-\frac{\partial\Psi^{e}}{\partial\mathbf{\epsilon} ^{e}}\right):\dot{\mathbf{\epsilon}}^{e}+\mathbf{\sigma}:\dot{\mathbf{\epsilon}}^{p}-R \dot{r}\geq 0 \tag{27}\]

Figure 13: Midbrain Section 2 response. (a) Training loss and number of active parameters over the training process, (b) Uniaxial tension data, the fit of the proposed approach, and fit of Ref. [31], (c) torsion data, the fit of the sparsified neural network approach, compared to fit of Ref. [31].
where \(R=\frac{\partial\Psi^{p}}{\partial r}\) is the thermodynamic force conjugate to \(r\). In order to guarantee the fulfillment of eq. (27) we set \(\mathbf{\sigma}=\frac{\partial\Psi^{e}}{\partial\mathbf{e}^{e}}\) and introduce the yield function \(f(\mathbf{\sigma},R)\) which is required to be a non-negative convex function of its arguments and zero-valued at the origin, i.e. \(f(\mathbf{0},0)=0\), see Ref. [69]. This function allows us to derive the evolution equations

\[\dot{\mathbf{\epsilon}}^{p}=\dot{\lambda}\frac{\partial f}{\partial\mathbf{\sigma}}, \quad\dot{r}=-\dot{\lambda}\frac{\partial f}{\partial\mathbf{R}} \tag{28}\]

with the consistency parameter \(\dot{\lambda}\). Many materials (such as metals or rocks) are characterized by pressure-independence and by linear elastic behavior in the elastic regime, i.e. where \(f(\mathbf{\sigma},R)<0\). The latter allows us to define the elastic component of the free energy function as

\[\Psi^{e}(\mathbf{\epsilon}^{e})=\frac{1}{2}\mathbf{\epsilon}^{e}:\mathbb{C}:\mathbf{ \epsilon}^{e} \tag{29}\]

where \(\mathbb{C}\) is a fourth-order material tangent specifying the (elastic) anisotropy. Under the assumption of isotropic yielding and due to the pressure independence of metals, the yield function can be rewritten as a function of two \(\pi\)-plane components

\[f(\frac{1}{R}\pi_{1},\frac{1}{R}\pi_{2}) \tag{30}\]

with

\[\begin{bmatrix}\pi_{1}\\ \pi_{2}\\ \pi_{3}\end{bmatrix}=\begin{bmatrix}\sqrt{\frac{2}{3}}&-\sqrt{\frac{1}{6}}&- \sqrt{\frac{1}{6}}\\ 0&\sqrt{\frac{1}{2}}&-\sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}\end{bmatrix}\begin{bmatrix} \sigma_{1}\\ \sigma_{2}\\ \sigma_{3}\end{bmatrix} \tag{31}\]

where \(\sigma_{i}\), \(i=1,2,3\) are the principal stresses and where \(R\) now acts as the ratio of a homothetic transformation of the initial yield function (\(R=1\)). For \(0<R<1\) the yield function is expanding and isotropic hardening can be represented. In the following, we will use the sparse regression framework to find representations for the initial yield function \(f(\pi_{1},\pi_{2})\) and for the hardening function \(R(r)\), respectively.

### Yield function from numerical data

We aim to represent an isotropic, pressure-independent initial yield function of the form \(f(\pi_{1},\pi_{2})\). Following Ref. [70] we can prove that for thermodynamic consistency the yield function is convex with regards to its inputs \(\pi_{1}\) and \(\pi_{2}\). We can therefore train a sparse input convex neural network (c.f. Section 2.1.1) and find an interpretable formulation for the yield functions. Since we focus on pressure-independent yield functions we introduce the deviatoric stress \(\mathbf{s}\) as well the two invariants \(J_{2}\) and \(J_{3}\), i.e.

\[\mathbf{s}=\mathbf{\sigma}-\frac{1}{3}\text{tr}(\mathbf{\sigma})\mathbf{I},\qquad J_{2}= \frac{1}{2}\text{tr}(\mathbf{s}^{2}),\qquad J_{3}=\frac{1}{3}\text{tr}(\mathbf{s}^{3}). \tag{32}\]

To make them easier to read the following representations are rounded to three decimal places. For all the examples the yield functions are trained on 30 training points, shown in Table A5.

### Drucker

The first example focuses on the well-known Drucker yield function [71] which involves both invariants of the Cauchy stress deviator and is here specified as

\[f_{D}=J_{2}^{3}+1.5\,J_{3}^{2}-(0.24)^{6}. \tag{33}\]

Figure 13(a) shows the training loss and the number of active parameters over the training process of the training run with the median final loss. We can see that the final model contains roughly \(10\) parameters at the end. The obtained, interpretable functional form reads

\[\begin{split}\hat{f}_{D}&=0.068\log\left(\left(e^{-2.122\pi_{1}+3.648\pi_{2}}+1\right)^{5.238}+1\right)+0.068\log\left(\left(e^{-2.121\pi_{1}-3.648\pi_{2}}+1\right)^{5.229}+1\right)\\ &+0.296\log\left(e^{4.922\pi_{1}}+1\right)-1.699.\end{split} \tag{34}\]

Using this form the approximated yield surface and the raw data are shown in Figure 13(b). We can see that the model fits the data proficiently well. Given the simplicity of the model we can furthermore, for example, use the series expansion around \(x=0\) of

\[\log\left(\frac{(e^{ax}+1)^{y}+1}{(e^{-ax}+1)^{y}+1}\right)\approx-\frac{a2^{y }yx}{2^{y}+1} \tag{35}\]
to show that

\[\begin{split} R&=\hat{f}_{D}(0,\pi_{2})-\hat{f}_{D}(0,- \pi_{2})\\ &=0.068\log\left(\frac{\left(e^{3.648\pi_{2}}+1\right)^{5.238}+1}{ \left(e^{-3.648\pi_{2}}+1\right)^{5.238}+1}\right)+0.068\log\left(\frac{\left( e^{-3.648\pi_{2}}+1\right)^{5.229}+1}{\left(e^{+3.648\pi_{2}}+1\right)^{5.229}+1} \right)\\ &\approx-0.068\frac{2.12225^{.238}5.238\pi_{2}}{2^{5.238}+1}+0.068 \frac{2.1212^{5.229}5.229\pi_{2}}{2^{5.229}+1}\\ &\approx-0.0018\pi_{2}\end{split} \tag{36}\]

the obtained yield function is roughly symmetric along \(\hat{f}_{D}(0,\pi_{2})\).

### Cazacu

The Drucker yield function is symmetric in the \(\pi_{1}\)-\(\pi_{2}\)-plane. To highlight that the model is also able to obtain an accurate model when tension-compression asymmetries are presented we consider the yield function suggested by Cazacu et. al. [72] which we specify to be

\[f_{C}=\left(|s_{1}|+0.5s_{1}\right)^{2}+\left(|s_{2}|+0.5s_{2}\right)^{2}+ \left(|s_{3}|+0.5s_{3}\right)^{2}-0.24. \tag{37}\]

We again obtain 30 data points on the yield limit of this model and fit a \(\mathcal{L}^{0}\)-regularized input convex neural network. Figure 14(a) shows the loss and parameter evolution. The functional form of the final model is given by

\[\begin{split}\hat{f}_{C}&=17.327\log\left(\left(1+ e^{-0.364\pi_{1}}\right)^{0.939}\left(1+e^{-0.313\pi_{2}}\right)^{1.852}+1 \right)\\ &+1.119\log\left(\left(e^{-7.341\pi_{1}+4.18\pi_{2}}+1\right)^{0. 95}e^{6.953\pi_{1}+3.93\pi_{2}}+1\right)-38.066.\end{split} \tag{38}\]

Figure 14(b) shows that the presented approach is also able to fit this asymmetric yield function with reasonable accuracy.

### Tresca yield function

In contrast to the other two yield functions, the last example deals with the non-smooth Tresca yield surface as proposed in Ref. [73] which we specify as

\[f_{T}=\max\left(|\sigma_{1}-\sigma_{2}|,|\sigma_{1}-\sigma_{3}|,|\sigma_{3}- \sigma_{2}|\right)-0.24. \tag{39}\]

Figure 14: Fit of a convex neural network to the Drucker yield function. (a) Training loss behavior and number of active (\(|w|>0\)) trainable parameters; (b)Data points (black dots), true curve (black dotted), and predicted curve (blue).
Figure (a)a again shows that both the number of parameters and the loss decrease over the course of the training process. The final functional representation is given by

\[\begin{split}\hat{f}_{T}&=0.021\log\left(\left(1+e^{- 116.395\pi_{2}}\right)^{2.706}+1\right)+0.019\log\left(\left(e^{-96.41\pi_{1}+55. 983\pi_{2}}+1\right)^{3.117}+1\right)\\ &+0.023\log\left(\left(e^{99.193\pi_{1}+57.03\pi_{2}}+1\right)^{2.46}+1\right)-1.127\end{split} \tag{40}\]

which is reasonably small to allow for some interpretation. Finally, the predicted yield limit \(\hat{f}_{T}=0\) is overlaid over the true data in Figure (b)b. We can see that the presented approach is also able to find interpretable functional forms when approximating non-smooth yield surfaces. We remark, that due to the smoothness of the activation function, the final yield function is necessarily also smooth. To exactly fit non-smooth yield surfaces other sharper activation functions like the rectified linear unit [74] could be employed.

Figure 16: Fit of a convex neural network to the Tresca yield function. (a) Training loss behavior and number of active (\(|w|>0\)) trainable parameters; (b)Data points (black dots), true curve (black dotted), and predicted curve (blue).

Figure 15: Fit of a convex neural network to the Cazacu yield function. (a) Training loss behavior and number of active (\(|w|>0\)) trainable parameters; (b)Data points (black dots), true curve (black dotted), and predicted curve (blue).
### Isotropic hardening law

For this application, we only have access to uniaxial monotonic loading tests and we aim to fit the isotropic hardening function \(R(r)\) as given in eq. (30). Given this limited data, we assume that the elastic response is isotropic, i.e.

\[\mathbb{C}=\frac{E\nu}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}+\frac{E}{2(1+\nu )}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right) \tag{41}\]

where the Young's modulus \(E\) and the Poisson's ratio \(\nu\) are adjustable parameters. We then aim to find the functional form \(R(r)\) depending on the internal hardening variable \(r\) which influences the scaling of a von Mises yield function [75] of the form

\[f(\frac{1}{R(r)}\,\pi_{1},\frac{1}{R(r)}\,\pi_{2})=\frac{3}{2}\sqrt{\frac{1}{R (r)}\,\pi_{1}^{2}+\frac{1}{R(r)}\,\pi_{2}^{2}}-\sigma_{y} \tag{42}\]

where the yield stress \(\sigma_{y}\) can be directly obtained from the data. If we assume isotropic hardening then the function \(R(r)\) is required to be positive and monotonically increasing. We can therefore employ a positive, monotonically increasing neural network as introduced in Section 2.1.2 with a sigmoid activation function to approximate the functional form of \(R(r)\). We remark that we require \(R(r=0.0)=1.0\) when no plastic yielding has occurred.

We will test the approach out on three different datasets which are summarized in Table 6.

#### 3.9.1 Experimental data - U71Mn rail steel

We start with a monotonic uniaxial loading dataset of U17Mn rail steel discussed in Ref. [76]. We first fit the elastic parameters and the yield stress to \(E=220\cdot 10^{3}\)MPa, \(\nu=0.3\) and \(\sigma_{y}=484.5\)MPa respectively. Using the network with the median loss after 10 runs, the training evolution and the reduction of the number of active parameters are shown in Figure 16(a). The obtained functional form for the isotropic hardening function reads

\[\hat{R}_{U71Mn}(r)=0.099+\frac{1.801}{1+e^{-194.688r}} \tag{43}\]

which is an easily interpretable function, i.e. we can see that it is always monotonically increasing since \(\exp(-r)\) is a monotonically decreasing function in \(r\). Figure 16(b) shows the true data, and the predicted response of the monotonic loading process with the fitted \(\hat{R}_{U71Mn}(r)\). The blue line indicates the range of the training data, i.e. \(R(r)\) was trained until roughly \(4\%\) strain. To highlight the extrapolation quality of the model the red line is the model prediction into unseen loading ranges. We can see that both the training and testing data were fitted proficiently well.

#### 3.9.2 Experimental data - SS316L stainless steel

Next, we look at data of a uniaxial loading test done on SS316L stainless steel performed by Ref. [76]. The fitted parameters of the elastic range read \(E=190\cdot 10^{3}\)MPa, \(\nu=0.35\) and the yield stress is given by \(\sigma_{y}=200\)MPa.

Figure 17: Fit of experimental monotonic loading curve with the hardening function represented by a sparse neural network - U71Mn rail steel. (a) Loss and number of active parameters over epochs (b) Uniaxial tension curve; Black dots represent the experimental data; Blue line indicates prediction in training data; Red line is extrapolation.
Figure (a)a depicts the evolution of the training loss and the parameters while Figure (b)b highlights the ability of the resulting model to fit the training data (up to the end of the blue solid line) and the ability to generalize well beyond the training region (solid red line). The interpretable functional form of the isotropic hardening function reads

\[\hat{R}_{SS316L}(r)=0.023+\frac{1.662}{1+1.071e^{-190.683r}}+\frac{0.362}{1+1.071 e^{-2200.640r}} \tag{44}\]

which, for example, allows us to easily see that there is no hardening when the internal variable is zero, i.e. \(\hat{R}_{SS316L}(0)\approx 0.023+0.802+0.175\approx 1\).

#### 3.9.3 Experimental data - 40Cr3MoV bainitic steel

Lastly, we look at monotonic loading data of 40Cr3MoV bainitic steel which was also published Ref. [76]. We set the material parameters as \(E=207\cdot 10^{3}\)MPa, \(\nu=0.3\) and \(\sigma_{y}=1000\)MPa. The loss and parameter evolution are shown in Figure (a)a. Figure (b)b depicts the interpolated and extrapolated predictions using the fitted model which is of the form

\[\hat{R}_{40Cr3MoV}(r)=-0.501+\frac{2.904}{1.669e^{-\frac{5.892}{1+1.720e^{-1 02.643r}}}\frac{0.740}{1+1.720e^{-667.227r}}+1}. \tag{45}\]

We can see that similar to the first two cases the accuracy of the fit is proficiently high.

## 4 Discussion and Conclusion

We have proposed to prune physics-augmented neural network-based constitutive models using a smoothed version of \(L^{0}\)-regularization to enable interpretable and trustworthy model discovery in a wide array of problems in mechanics. The network is trained by simultaneously fitting the training data and penalizing the number of active parameters. On a variety of applications including synthetic and experimental data, we have shown that we are able to obtain accurate, yet interpretable constitutive models for compressible and incompressible hyperelasticity, yield functions, and isotropic hardening functions. The presented approach seems to be highly flexible and has the potential to overcome the restrictions of functional form selection, towards automation of constitutive modeling. Uniquely, i) enforcing physical constraints enables generalization/extrapolation as well as training with limited and low data, ii) utilization of neural networks enables high expressiveness and eliminates the development of specific model form libraries, and iii) pruning, leads to interpretable discovery and also enhances generalization/extrapolation.

In the next steps, we will use this approach to obtain functional forms for other constitutive models such as representations for kinematic hardening or viscoelasticity. We have furthermore (purposefully) only used neural networks

Figure 18: Fit of experimental monotonic loading curve with the hardening function represented by a sparse neural network - SS316L stainless steel. (a) Loss and number of active parameters over epochs (b) Uniaxial tension curve; Black dots represent the experimental data; Blue line indicates prediction in training data; Red line is extrapolation.
with single nonlinear activation functions. Future work could center around finding representations using a different activation function at each neuron as a means to enhance the expressivity of our representations to parsimoniously model even more complex material responses.

## Acknowledgments

JF and NB gratefully acknowledge support by the Air Force Office of Scientific Research under award number FA9550-22-1-0075.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.
