A learning-based optimal uncertainty quantification method and its application to ballistic impact problems

We are concerned with a system whose performance is described by a known response function (1)$Y=F\left(X\right),$from one compact measurable space $\mathcal{X}\subseteq {\mathbb{R}}^m$ of inputs to a second measurable space $\mathcal{Y}\subseteq {\mathbb{R}}^n$ of outputs. $X\equiv (x_1,\dots ,x_m)$ are $m$ real-valued random variables, expressing imperfectly known or uncertain properties of the system and $Y\equiv (y_1,\dots ,y_n)$ are $n$ real-valued random variables, or performance measures. The input variables are generated randomly according to an unknown random variable $X$ with values in $\mathcal{X}$ according to a law $\mathbb{P}\in M\left(\mathcal{X}\right)$, where $M\left(\mathcal{X}\right)$ is the set of all probability measures supported on $\mathcal{X}$. The design specifications require that $Y$ remains that $Y\in {\mathcal{Y}}_{\text{a}}$ for some admissible set ${\mathcal{Y}}_{\text{a}}\subseteq \mathcal{Y}$. Thus the design fails if $Y\in {\mathcal{Y}}_{\text{c}}$ where ${\mathcal{Y}}_{\text{c}}=\mathcal{Y}\setminus {\mathcal{Y}}_{\text{a}}$ is the inadmissible set. Ideally, we would like the support of the probability measure associated to $Y$ to be contained within ${\mathcal{Y}}_{\text{a}}$, i.e., (2)$\mathbb{P}[Y\in {\mathcal{Y}}_{\text{a}}]=1.$Systems satisfying this condition can be certified with complete certainty. However, this absolute guarantee of safe performance may be unattainable, e.g., if $\mathbb{P}$ lacks compact support, or is prohibitively expensive. In these cases, we may relax the condition of certification. Let $ϵ\in [0,1]$ denote the greatest acceptable probability of failure (PoF). Then we say that the system is safe if (3)$\mathbb{P}[Y\in {\mathcal{Y}}_{\text{c}}]\le ϵ,$and the system is unsafe if (4)$\mathbb{P}[Y\in {\mathcal{Y}}_{\text{c}}]>ϵ.$

However, due to the lack of information, the exact probability measure $\mathbb{P}$ may be unknown. To this end, we define a subset (5)$\mathcal{A}\subseteq \left\{\mu \mid \mu \in M\left(\mathcal{X}\right)\right\}$that encodes all the information that we have about the probability measure of the random variables $\mathbb{P}$. This information may come from experimental data, lower-level simulations or expert opinions. Notably, $\mathbb{P}\in \mathcal{A}$ and some admissible scenarios in $\mathcal{A}$ may be safe (i.e., $\mu [Y\in {\mathcal{Y}}_{\text{c}}]\le ϵ$), whereas other admissible scenarios may be unsafe (i.e., $\mu [Y\in {\mathcal{Y}}_{\text{c}}]>ϵ$). Now observe that, given such an information/assumptions set $\mathcal{A}$, there exist upper and lower bounds on $\mathbb{P}[Y\in {\mathcal{Y}}_{\text{c}}]$ corresponding to the scenarios compatible with assumptions, i.e., the values $U\left(\mathcal{A}\right)$ and $L\left(\mathcal{A}\right)$ of the optimization problems (6a)$U\left(\mathcal{A}\right)≔\underset{\mu \in \mathcal{A}}{sup}\mu [Y\in {\mathcal{Y}}_{\text{c}}],$ (6b)$L\left(\mathcal{A}\right)≔\underset{\mu \in \mathcal{A}}{inf}\mu [Y\in {\mathcal{Y}}_{\text{c}}].$ Assume that $U\left(\mathcal{A}\right)$ and $L\left(\mathcal{A}\right)$ can indeed be computed on demand. Now, since $\mathbb{P}\in \mathcal{A}$, it follows that (7)$L\left(\mathcal{A}\right)\le \mathbb{P}[Y\in {\mathcal{Y}}_{\text{c}}]\le U\left(\mathcal{A}\right).$The upper bound $U\left(\mathcal{A}\right)$ is optimal in the sense that for all $\mu \in \mathcal{A}$, we have $\mu [Y\in {\mathcal{Y}}_{\text{c}}]\le U\left(\mathcal{A}\right)$, and if $U^{\prime }<U\left(\mathcal{A}\right)$, there exists $\mu \in \mathcal{A}$ such that $U^{\prime }<\mu [Y\in {\mathcal{Y}}_{\text{c}}]\le U\left(\mathcal{A}\right)$. Similar conclusions apply for the lower bound $L\left(\mathcal{A}\right)$.

The bounds $U\left(\mathcal{A}\right)$ and $L\left(\mathcal{A}\right)$ defined in Eq. (6) can be used to construct a solution to the certification problem. Provided that the information set $\mathcal{A}$ is valid (in the sense that $\mathbb{P}\in \mathcal{A}$), then if $U\left(\mathcal{A}\right)\le ϵ$, then the system is provably safe; if $ϵ<L\left(\mathcal{A}\right)$, then the system is provably unsafe; and if $L\left(\mathcal{A}\right)\le ϵ<U\left(\mathcal{A}\right)$, then the safety of the system cannot be decided due to lack of information. The corresponding certification process and its optimality are illustrated in Fig. 1. Evidently, the tighter the bounds $U\left(\mathcal{A}\right)$ and $L\left(\mathcal{A}\right)$ the more economical the certification. However, increasing tightness comes at increasing computational expense, which sets forth a fundamental trade-off between economy of certification and computability.

Our goal is to compute the optimal probability $U\left(\mathcal{A}\right)$ and $L\left(\mathcal{A}\right)$, as defined in Eq. (6). In general, the OUQ problem is an infinite dimensional optimization problem and is computationally intractable. To this end, we follow Winkler et al., 1979, Winkler, 1988 and show that for the constraint set $\mathcal{A}$ that is of practical interest in uncertainty quantification, the infinite-dimensional optimization problem can be reduced to a finite dimensional optimization problem. In this section, we list the main theorem from Winkler (1988), as well as its application to the OUQ problem. We refer the readers to Owhadi et al. (2013) for further references.

For the sake of brevity, we shall restrict our attention to $U\left(\mathcal{A}\right)$ as defined in Eq. (6), and the treatment for $L\left(\mathcal{A}\right)$ follows a similar manner. First note the set of all probability measures supported on $\mathcal{X}$, $M\left(\mathcal{X}\right)$, is convex. Indeed, given two arbitrary measures ${\mu }_1,{\mu }_2\in M\left(\mathcal{X}\right)$, the convex combination ${\mu }_3=\alpha {\mu }_1+(1-\alpha ){\mu }_2,\alpha \in (0,1)$ is also a probability measure in $M\left(\mathcal{X}\right)$. Furthermore, it is known that the extremal points of $M\left(\mathcal{X}\right)$ are Dirac measures of the form (Rudin, 1991) (8)${\delta }_z\left(\mathcal{Z}\right)=\begin{array}{cc}1\hfill & \text{if}z\in \mathcal{Z},\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$

Theorem 1 provides a path to reduce the infinite dimensional optimization problem defined in Eq. (6) to a finite dimensional optimization problem on the extremal points $\mathbf{ex}\mathcal{A}$.

To see this, first recall that a Choquet type integral representation formula can be derived for probability measures (Winkler, 1988). Specifically, for every probability measure in the constrained set $\mu \in \mathcal{A}$ defined in Theorem 1, there is a probability measure $p$ supported on $\mathbf{ex}\mathcal{A}$, such that (11)$\mu (Y\in {\mathcal{Y}}_c)={\int }_{\nu \in \mathbf{ex}\mathcal{A}}\nu (Y\in {\mathcal{Y}}_c)dp\left(\nu \right).$Subsequently, (12)$sup\{\mu (Y\in {\mathcal{Y}}_c):\mu \in \mathcal{A}\}=sup\left\{{\int }_{\mathbf{ex}\mathcal{A}}\nu \left(B\right)dp\left(\nu \right)\right\}=sup\{\nu (Y\in {\mathcal{Y}}_c):\nu \in \mathbf{ex}\mathcal{A}\}.$ Therefore, the calculation of $U\left(\mathcal{A}\right)$ can be reduced into an optimization problem of finding the supremum of extremal points in the constraint set $\mathcal{A}$, where (13)$U\left(\mathcal{A}\right)=\underset{\genfrac{}{}{0ex}{}{X_i\in \mathcal{X},i=1,\dots ,K+1}{t_i\in [0,1],i=1,\dots ,K+1}}{max}\sum _{i=1}^{K+1}{t}_i{\delta }_{F\left(X_i\right)}\left({\mathcal{Y}}_c\right),$$\text{subject to}\sum _{i=1}^{K+1}{t}_i=1,\sum _{i=1}^{K+1}{t}_i\cdot f_j\left(X_i\right)=c_j,j=1,\dots ,K.$

In addition, we note for the special case where the admissible set contains no constraint, the optimization problem reduced into optimizing a single Dirac measure, (14)$U\left(\mathcal{A}\right)=\underset{X\in \mathcal{X}}{max}{\delta }_{F\left(X\right)}\left({\mathcal{Y}}_c\right),$which corresponds to finding the worst case scenario. Furthermore, Theorem 1 provides a way of constraining the input/prior distribution. As an example, the mean constraint can be imposed by setting $f=X$, while the higher-order moment constraints can be imposed by setting $f=X^k$, where $k$ is the order of moment. We shall demonstrate this further in the example studies in Section 3.

The implementation of the optimization problem defined in Eq. (13) requires repetitive evaluation of the performance indicator: (15)$\mathcal{D}:X\mapsto {\delta }_{F\left(X\right)}\left({\mathcal{Y}}_c\right),$which is a composition of the system response function $Y=F\left(X\right)$ and the Dirac measure ${\delta }_Y\left({\mathcal{Y}}_c\right)$. Direct evaluation of the system’s response function can become very expensive for high-dimensional, non-linear systems. Our idea is to learn the performance indicator from data using a deep neural network and by utilizing data generated by solutions of the system’s response function over various inputs sampled from an appropriate probability measure defined in the space of input $\mathcal{X}$. To do so, we observe that the output of the performance indicator $\mathcal{D}$ can only take the values from the set $\{0,1\}$, hence the problem closely resembles the classification problem in machine learning (Osisanwo et al., 2017). To that end, we consider a neural network approximation of ${\mathcal{D}}_{NN}:\mathcal{X}\times {\mathbb{R}}^{n_{\theta }}\mapsto \{0,1\}$ with parameterization $\theta \in {\mathbb{R}}^{n_{\theta }}$, such that (16)${\mathcal{D}}_{NN}(X;\theta )\approx \mathcal{D}\left(X\right).$The proposed approach can then be finished in three steps:

We proceed with a series of comments on the learning-based OUQ approach as well as its alternatives.

We assume that the constitutive behavior of the plate is characterized by an appropriately calibrated Johnson-Cook plasticity model(Johnson, 1983), (19)$\sigma \left({ϵ}_p,{\dot{ϵ}}_p,T\right)=\left[A+B{ϵ}_p^n\right]\left[1+Cln{\dot{ϵ}}_p^{\ast }\right]\left[1-{T^{\ast }}^m\right],$where $\sigma$ is the true Mises stress, ${ϵ}_p$ is the equivalent plastic strain, ${\dot{ϵ}}_p$ is the plastic strain rate, and $T$ is the temperature. The normalized plastic strain rate ${\dot{ϵ}}_p^{\ast }$ and temperature $T^{\ast }$ are defined as (20)${\dot{ϵ}}_p^{\ast }≔\frac{{\dot{ϵ}}_p}{{\dot{ϵ}}_{p0}},$and (21)$T^{\ast }≔\frac{T-T_0}{T_m-T_0},$respectively, where ${\dot{ϵ}}_{p0}$ is a reference strain rate, $T_0$ is a reference temperature and $T_m$ is the melting temperature. The model parameters are: $A$, the yield stress; $B$, the strain-hardening modulus; $n$, the strain-hardening exponent; $C$, the strengthening coefficient of strain rate; and $m$, the thermal-softening exponent.

We choose our system’s performance measure $Y$, to be the maximum backface deflection of the plate after the impact $Y=y_r\in \mathbb{R}$. We set a maximum backface deflection threshold, $Y_T\in \mathbb{R}$, and consider the system safe if the $y_r\in {\mathcal{Y}}_a=[0,Y_T)$, and unsafe if $y_r\in {\mathcal{Y}}_c=[Y_T,+\infty )$. Furthermore, we regard the set $\{A,B,n,C,m\}$ of Johnson-Cook parameters as the main source of uncertainty in the analysis. The bounds of these input uncertain parameters are tabulated in Table 1, which are determined by experimental characterization with $95$% confidence intervals (Hasenpouth, 2010). In what follows, we normalize the Johnson-Cook parameters into the range $[0,1]$ using the min–max feature scaling and the corresponding fixed bounds in Table 1. As a result, the random inputs in the calculations are the normalized Johnson-Cook parameters represented by $X\equiv \{\bar{A},\bar{B},\bar{n},\bar{C},\bar{m}\}\in {[0,1]}^5$. We therefore aim to estimate the probability of failure (PoF) $\mathbb{P}[Y\in {\mathcal{Y}}_c]$, from uncertainties in inputs $X$.

We use finite element software LS-DYNA (Hallquist et al., 2007) to solve the impact problem. A schematic of the finite-element model is shown in Fig. 2. The diameter of the projectile is 1.12 cm, and the size of the plate is $10\times 10\times 0.35$ cm. The attack velocity is $200\mathrm{m/s}$ with normal impact. The backface nodes of the target near the edges are fully constrained to prevent displacement in all directions. The projectile is resolved using $864$ elements, while the number of elements for the plate is 70,000. All the elements are linear hex, single point integration with careful hourglass control. The time-step size is adaptive and determined by the critical size of elements, with all simulations running for $500.0\mathrm{\mu }\text{s}$ before termination. This simulation duration is sufficiently long to allow for the rebound and separation of the projectile from the plate in all the calculations. The calculations are adiabatic with the initial temperature set at room temperature. The equation-of-state, which controls the volumetric response of the material, is assumed to be of the Gruneisen type. For simplicity, the projectile is assumed to be rigid and uncertainty-free. All other material parameters are fixed and listed in Table 2.

We follow the learning-based OUQ method as detailed in Section 2.3. A total of 2000 input–output pairs $\{X,Y\}$ are generated using the Latin Hypercube Sampling method with $L=2000$ and $N_L=1$. We then use this data to compute the performance indicator $\mathcal{D}\left(X\right)$ for a variety of $Y_T\in [0.9,1.3]$ cm. We subsequently train a series of neural networks with the architecture described in Section 2.3. For each $Y_T$, we use a total of $1500$ samples to train and the remaining $500$ to test the learned indicator. In all cases, neural network consists of $4$ intermediate layers with $200$ nodes per layer and are trained using the ADAM (Kingma and Ba, 2014) method.

Fig. 3 shows the results of a typical learned performance indicator with threshold $Y_T=1.03$ cm. Both training and testing error is included in Fig. 3(a) with minimum training and testing error to be 0.005 and 0.0287 respectively. A set of $50$ testing samples were randomly selected with the input $X=\{\hat{A},\hat{B},\hat{n},\hat{C},\hat{m}\}$ plotted in Fig. 3(b). The resultant output ${\delta }_{F\left(X\right)}\left({\mathcal{Y}}_c\right)$ is plotted in Fig. 3(c). It is evident from Fig. 3(c) that the neural network prediction (marked in square) agrees very well with the true solution (marked in circle) regarding the cases of both $0$ and $1$.

We are now in a position to demonstrate the power of our method by considering a series of case studies with varying levels of knowledge on input probabilities. We recall that such knowledge is applied as constraints in the optimization problems, Eq. (13). In this section, we will study multiple cases including mean and higher-order moment constraints, and complete and partial moment constraints. The corresponding constraints are also expressed as mathematical equations in our OUQ framework. In all cases, we have tested the penalty coefficients ${\lambda }_i,i=1,\dots ,K$, over multiple orders of magnitude ranging from $10$ to $10^5$. We finally fix all ${\lambda }_i$ at $10^3$ which provides an excellent compromise between computational speed and constraint satisfaction. The computation time required to solve the optimization problems is on the order of seconds.

The dependence of the upper and lower PoF bounds on the means of $\bar{A}$ and $\bar{B}$ is shown in Fig. 7. These bounds are determined by solving the OUQ optimization at different mean values of $\bar{A}$ or $\bar{B}$ while constraining the other four normalized Johnson-Cook parameters in the corresponding ranges with a fixed mean of 0.5. The threshold of the maximum backface deflection is also fixed at $1.03$ cm. As expected, for both $\bar{A}$ and $\bar{B}$, the PoF bounds decrease with the increase of means of $\bar{A}$ and $\bar{B}$, since the plate is more likely to be stronger and stiffer. It is also notable in Fig. 7(a) that when the mean of $\bar{A}$ is less than 0.09 (i.e., $\mathbb{E}\left(A\right)=205\mathrm{MPa}$), both the upper and lower bounds on the PoF converge to $1$. This means that the system is unsafe with 100% confidence regardless of the values of other material properties as long as they satisfy the constraints. Therefore, under this scenario the impact performance of the plate is governed by the yield stress of AZ31B Mg alloys. Moreover, each subfigure of Fig. 7 is divided by the lower and upper bounds into three regions labeled as “Decertified”, “Cannot decide” and “Certified”. As a result, for a given tolerance $ϵ$ on the PoF, these regions provide reference ranges of the design parameters to achieve the rigorous certification criterion illustrated in Fig. 1. That is, the system is provable safe if the design parameter and the PoF tolerance are located in the region “Certified”, is provable unsafe if in the region “Decertified”, and the safety of the system cannot be decided if in the region “Cannot decide”.