On Taylor Series Expansion for Statistical Moments of Functions of Correlated Random Variables ^†

Abstract

The paper is focused on Taylor series expansion for statistical analysis of functions of random variables with special attention to correlated input random variables. It is shown that the standard approach leads to significant deviations in estimated variance of non-linear functions. Moreover, input random variables are often correlated in industrial applications; thus, it is crucial to obtain accurate estimations of partial derivatives by a numerical differencing scheme. Therefore, a novel methodology for construction of Taylor series expansion of increasing complexity of differencing schemes is proposed and applied on several analytical examples. The methodology is adapted for engineering applications by proposed asymmetric difference quotients in combination with a specific step-size parameter. It is shown that proposed differencing schemes are suitable for functions of correlated random variables. Finally, the accuracy, efficiency, and limitations of the proposed methodology are discussed.

1. Introduction

Mathematical modeling in civil engineering is often represented by the finite element method (FEM). Although FEM is an accurate and efficient technique, it is still highly time-consuming, particularly in the case of non-linear FEM including geometrical and material non-linearity. Therefore, from a practical point of view, it is necessary to decrease the number of FEM calculations as much as possible while satisfying the given safety requirements of the analyzed structure. A solution can be represented by a semi-probabilistic approach widely accepted in the engineering field [1] and implemented into the national codes such as Eurocode [2]. Such approach is able to greatly reduce the number of necessary calculations for the design and an assessment of structures. The basic reliability concept is given as $Z=R-E$, where Z is a safety margin, which is defined as the difference between the structural resistance R and the load effect E. The task of reliability analysis is the estimation of failure probability $p_f=P\left(Z<0\right)$, which might be highly computationally demanding. According to the semi-probabilistic approach, the resistance of a structure R is separated, and the design value $R_d$ satisfying given safety requirements is evaluated instead of calculating the failure probability. Such approach directly leads to the design value of resistance, which is obtained by the traditional Partial Safety Factor (PSF) approach, and thus can be easily used for a design and an assessment of structures. The PSF method is based on a simple assumption, that a calculation with design values of input random variables leads to the design value of resistance $R_d=r\left({\mathbf{x}}_{\mathbf{d}}\right)$, where design values of input random variables ${\mathbf{x}}_{\mathbf{d}}$ are derived under several simplifications, such as a linearization of a limit state function. In consequence, PSF works well for standard linear calculations, but there may be a significant error for a non-linear analysis, which is far more popular nowadays. Therefore, it is necessary to develop new methods in compliance with the semi-probabilistic approach applicable for non-linear analysis. The semi-probabilistic approach is briefly presented in the following paragraph.

It is assumed that R and E are independent, and separated R is lognormally distributed; thus, the design value of resistance $R_d$ is defined as

$$R_d={\mu }_R·exp(-{\alpha }_R\beta v_R),$$

(1)

where $v_R$ is the coefficient of variation (CoV) of resistance, and ${\alpha }_R$ represents the sensitivity factor associated with R derived from the First Order Reliability Method (FORM) [1,3]. FORM is commonly applied to linearization of limit state function at the most probable failure point by Taylor series expansion. FORM assumes the uncorrelated standardized Gaussian space $\xi$; thus, all variables must be transformed by Rosenblatt transformation [4] from the original space. The coordinates of the most probable failure point, also called the design point, are thereafter described by the shortest distance $\beta$ to the origin of the $\xi$ space, direction cosines ${\alpha }_R$ associated with resistance and ${\alpha }_E$ associated with load. The shortest distance $\beta$ is defined as the Hasofer–Lind reliability index, and its minimal value is given for various conditions in normative documents, in order to achieve the target safety of structures.

For industrial applications, FORM is simplified by the statistical estimation of fixed value ${\alpha }_R=0.8$. Therefore, to determine the design value by a semi-probabilistic approach, it is crucial to correctly estimate the mean value and variance of structural resistance R, which can be seen as a function of multiple random variables. This task may be challenging due to the fact that input random variables can generally be non-Gaussian and correlated. There have been several methods proposed in last two decades to estimate the variation coefficient of R (ECoV methods) [5,6,7,8,9,10]; however, mathematical background and limitations of these methods are often missing, and there is no solution for correlated random variables, which is common for material characteristics.

The only general approach to estimate statistical moments is pseudo-random sampling by a Monte Carlo type algorithm such as Crude Monte Carlo or Latin Hypercube Sampling [11,12] employed in numerical examples as a reference solution. However, it is necessary to perform a high number of simulations of the original mathematical model, which is not feasible in industrial applications due to the enormous computational burden. On the other hand, it is possible to assume several simplifications and create an approximation of the original mathematical model of R.

The approximating function is called a surrogate model, or a metamodel, and it is a topic of great interest among researchers from various research fields. The Polynomial Chaos Expansion (PCE) is often used for uncertainty quantification [13,14]. The Gaussian process or Krigging has recently received significant attention in reliability analysis of systems with very low failure probabilities [15], and artificial neural networks are often utilized for reliability-based optimization [16]. Although PCE, Krigging, and ANN represent very powerful and efficient approaches with many advantages, these advanced techniques require deep knowledge of theoretical background, and it is necessary to use developed algorithms with great caution.

Another well-known approximation of functions is Taylor series expansion (TSE), which was also used for derivation of PSF and FORM as described, for example, in [17]. Although TSE is often used to estimate statistical moments of functions of random variables by mathematicians, it has not yet been well investigated in the context of non-linear FEM in civil engineering in order to adapt and directly use TSE for structural reliability and semi-probabilistic approaches. For industrial applications, it is crucial that the proposed methods are easy to implement and easy to use with the same level of knowledge about the mathematical model as in the case of PSF. Therefore, it makes perfect sense to generalize TSE, which is already utilized for the derivation of PSF implemented in codes, to directly use in combination with FEM and semi-probabilistic ECoV approach. Therefore, the ECoV method based on TSE adapted for civil engineers is discussed, and several modifications of this approach are proposed in the next section. Moreover, the whole methodology of increasing complexity and accuracy of TSE suitable for industrial applications is proposed in this paper.

2. Taylor Series Expansion

An original mathematical model is often highly time-consuming, and it is necessary to create an approximation—a simplified function in explicit form. Although there are several advanced types of surrogate models, it is still common to use the traditional approach, called Taylor series expansion, which can be truncated to arbitrary order and used with various differencing schemes. Although such adaptivity makes TSE a powerful technique, there are severe problems for practical computations in the case of non-linear functions with complex stochastic models containing a dependence structure. In the following paragraphs, let us assume an original mathematical model in form of software algorithm (e.g., FEM); thus, the derivatives must be calculated numerically.

Let ($\mathsf{\Omega },\mathcal{F},\mathcal{P}$) be a probability space, where $\mathsf{\Omega }$ is an event space, $\mathcal{F}$ is Borel $\sigma$-algebra on $\mathsf{\Omega }$, and $\mathcal{P}$ is a probability measure $\mathcal{P}:\mathcal{F}\to \left[0,1\right]$. Let us assume a random vector $\mathbf{X}={\left(X_1,X_2,\dots X_n\right)}^T$ consisting of random variables $X\left(\omega \right),\omega \in \mathsf{\Omega }$ with existing mean values ${\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}$ and a mathematical model of this input random vector $r\left(\mathbf{X}\right)$. The response of the mathematical model is thereafter a random variable R described by a specific probability distribution and statistical moments. Further, let us assume the mathematical model $r\left(\mathbf{X}\right)$ to be infinitely differentiable in some open interval around the vector of mean values ${\mu }_{\mathbf{X}}={\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}$. Under this assumption, it is possible to expand the original model to the infinite Taylor series according to Taylor’s theorem:

$$\begin{array}{c}r\left(\mathbf{X}\right)=r\left({\mathit{\mu }}_{\mathit{X}}\right)+\nabla r\left({\mathit{\mu }}_{\mathit{X}}\right)·(\mathbf{X}-{\mathit{\mu }}_{\mathit{X}})+\frac{1}{2}(\mathbf{X}-{\mathit{\mu }}_{\mathit{X}})·\nabla \nabla r\left({\mathit{\mu }}_{\mathit{X}}\right)·(\mathbf{X}-{\mathit{\mu }}_{\mathit{X}})+···=\hfill \\ \hfill =r\left({\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}\right)+\sum _{i=1}^n\frac{\partial r\left(X\right)}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\frac{{\partial }^{2}r\left(X\right)}{\partial X_i{X}_j}\left(X_i-{\mu }_{X_i}\right)\left(X_j-{\mu }_{X_j}\right)+\dots \end{array}$$

(2)

where the derivatives are evaluated at ${\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}$. Note that TSE consists of a constant term, linear term, quadratic term, etc. For a practical computation, it is crucial to reduce Taylor series to a finite number of terms and to obtain derivatives by numerical differentiation. There are many possible differencing schemes, which are more or less suitable for specific applications. One of the possible formulas for numerical derivation was proposed by Schlune et al. [9], especially for civil engineers, where derivatives are approximated by the asymmetric difference quotient as follows:

$$\frac{\partial r\left(X\right)}{\partial X_i}=\frac{R_{X_m}-R_{X_{i\Delta }}}{{\Delta }_{X_i}}.$$

(3)

where the response of mathematical model $R_{X_m}$ is a calculation with mean values of X, and $R_{X_{i\Delta }}$ is the result of the model using reduced mean values of the i-th input random variables by ${\Delta }_{X_i}$. This differencing scheme is adapted for a structural design and an assessment by the step-size parameter $c=\left({\alpha }_R\beta \right)/\sqrt{2}$, and $X_{i\Delta }$ corresponds to quantile $F_i^{-1}\left(\mathsf{\Phi }(-c)\right)$, where $F_i^{-1}$ is an inverse cumulative distribution function of the i-th variable, and $\mathsf{\Phi }$ is the cumulative distribution function of standardized Gaussian distribution. For the sake of clarity, the difference is calculated as ${\Delta }_{X_i}=X_{im}-X_{i\Delta }$. Note that the step-size parameter is a function of the reliability index; thus, it is in compliance with the philosophy of a semi-probabilistic approach implemented in civil engineering codes [1]. Following this idea, additional asymmetric differencing schemes adapted for civil engineering used in combination with TSE of the first and the second order are proposed in the following subsections.

2.1. Linear Terms of Taylor Series Expansion

In engineering applications, it is common to assume only linear terms of TSE and independent input random variables. Since a semi-probabilistic approach is focused on practical applications, the significant advantage of TSE reduced to linear terms is the possibility of analyzing expressions for an expected value and a variance, see e.g., [17].

Theorem 1.

If an original mathematical function $r:{\mathbb{R}}^n\to \mathbb{R}$ of n independent random variables described by mean value ${\mu }_{X_i}$ and variance ${\sigma }_{X_i}^2$ is approximated by Taylor series expansion reduced to linear terms, the first two statistical moments of the response $R_T$ of linear Taylor approximation are analytically obtained as follows:

$${\mathbb{E}}_{R_T}\approx r\left({\mathit{\mu }}_{\mathit{X}}\right)$$

(4)

$$Va{r}_{R_T}\approx \sum _{i=1}^n{\left(\frac{\partial r\left(X\right)}{\partial X_i}\right)}^2{\sigma }_{X_i}^2$$

(5)

Proof of Theorem 1.

For the sake of clarity, the estimations of expected value $E_{R_T}$ and variance $Va{r}_{R_T}$ for the function of n independent random variables are as follows:

$${\mathbb{E}}_{R_T}\approx \mathbb{E}\left[r\left({\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}\right)\right]+\sum _{i=1}^n\mathbb{E}\left[\frac{\partial r\left(X\right)}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)\right]\approx r\left({\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}\right)$$

(6)

and

$$\begin{array}{c}Va{r}_{R_T}\approx Var\left[r\left({\mu }_{X_1},{\mu }_{X_2},\dots ,{\mu }_{X_n}\right)+\sum _{i=1}^n\frac{\partial r\left(X\right)}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)\right]=\sum _{i=1}^{n}Var\left[\frac{\partial r\left(X\right)}{\partial X_i}\left(X_i-{\mu }_{X_i}\right)\right]+\hfill \\ \hfill +\underset{i\ne j}{\sum _{i,j=1,\dots ,N}}Cov\left[\frac{\partial r\left(X\right)}{\partial X_i}\left(X_i-{\mu }_{X_i}\right),\frac{\partial r\left(X\right)}{\partial X_j}\left(X_j-{\mu }_{X_j}\right)\right]=\sum _{i=1}^n{\left(\frac{\partial r\left(X\right)}{\partial X_i}\right)}^2{\sigma }_{X_i}^2\end{array}$$

(7)

where the final equation arises from the definition of variance $Var\left(X\right)={\sigma }_X^2=\mathbb{E}\left[{\left(X-\mu \right)}^2\right]$ and property of variance $Var(c{X}_i+d{X}_j)=c^{2}Var\left(X_i\right)+d^{2}Var\left(X_j\right)+2cdCov(X_i,X_j)$. Moreover, for independent variables, the covariance between variables is equal to zero, and thus the formula is reduced.  □

As can be seen from the proof above, there is a strict assumption of uncorrelated random variables for Equation (5). However, it is necessary to assume correlated random variables in some practical examples solved by FEM to represent realistic behaviors of structures. An extension of the method for dependent random variables can be obtained from the proof above using first-order Taylor series expansion assuming correlation among random variables represented by the correlation coefficient $\rho$ in analytical form as

$$Va{r}_{R_T}\approx \sum _{i=1}^n{\left(\frac{\partial r\left(X\right)}{\partial X_i}\right)}^2{\sigma }_{X_i}^2+\underset{i\ne j}{\sum _{i,j=1,\dots ,n}}{\rho }_{i,j}{\sigma }_{X_i}{\sigma }_{X_j}\frac{\partial r\left(X\right)}{\partial X_i}\frac{\partial r\left(X\right)}{\partial X_j}.$$

(8)

However, higher terms of TSE or more accurate approximation of derivatives should be considered for the correct estimation of variance in the case of dependent input random variables and non-linear functions. Otherwise, the correlation term may lead to significant inaccuracy of the resulting variance. We propose the second-order backward asymmetric differencing according to Equation (9), which is adapted for structural design utilizing the parameter $c=\left({\alpha }_R\beta \right)/\sqrt{2}$ analogously to Equation (3) proposed by Schlune et al. The middle additional term $R_{Xi\frac{\Delta }{2}}$ is obtained by an evaluation of the original mathematical model with reduced i-th variable $X_{i\frac{\Delta }{2}}=X_{im}-{\Delta }_{X_i}/2$. Note that the proposed approach needs $2n+1$ evaluations of the original model, while the scheme proposed by Schlune needs $n+1$ simulations. In practice, an analyst could use the derivative scheme according to Schlune and further compute additional n simulations in order to obtain $R_{Xi\frac{\Delta }{2}}$ and more accurate results.

$$\frac{\partial r\left(X\right)}{\partial X_i}=\frac{3R_{X_m}-4R_{Xi\frac{\Delta }{2}}+R_{X_{i\Delta }}}{{\Delta }_{X_i}}.$$

(9)

2.2. Higher-Order Taylor Series Expansion

If higher terms of TSE are considered, it is inefficient to derive analytical formulas for statistical moments [18], and thus mean and variance should be calculated numerically by simulation techniques directly from Equation (2) truncated to quadratic terms. Moreover, additional higher-order derivatives must be evaluated, which might not be feasible in computationally demanding practical examples. Therefore, linear TSE is preferred for practical computations. However, for specific cases with significant interaction of input variables, one may use second-order TSE for the estimation of coefficient of variation. In this case, it is necessary to compute all second-order partial derivatives. For numerical calculations of $\frac{{\partial }^{2}r\left(X\right)}{\partial X_i^2}$, it is possible to use the already defined simulations $R_{X_m},R_{X_{i\Delta }},R_{Xi\frac{\Delta }{2}}$ in a standard asymmetric backward differencing scheme:

$$\frac{{\partial }^{2}r\left(X\right)}{\partial X_i\partial X_i}=\frac{R_{X_m}-2R_{Xi\frac{\Delta }{2}}+R_{X_{i\Delta }}}{{\Delta }_{X_i}^2}$$

(10)

The only additional computations of the original mathematical model needed are for mixed partial derivatives $\frac{{\partial }^{2}r\left(X\right)}{\partial X_i\partial X_j}$. Note that it is necessary to perform additional $\left(\genfrac{}{}{0pt}{}{n}{2}\right)$ simulations in order to obtain all the mixed partial derivatives. In total, it is necessary to calculate $2n+\left(\genfrac{}{}{0pt}{}{n}{2}\right)+1$ simulations for second-order TSE using the proposed asymmetric differencing schemes.

Theorem 2.

Mixed partial derivatives can be approximated by the simple backward finite differencing as

$$\frac{{\partial }^{2}r\left(X\right)}{\partial X_i\partial X_j}=\frac{R_{X_m}-R_{X_{i\Delta }}-R_{X_{j\Delta }}+R_{X_{i\Delta }{X}_{j\Delta }}}{{\Delta }_{X_i}{\Delta }_{X_j}},$$

(11)

where $R_{X_{i\Delta }{X}_{j\Delta }}$ represents the response of a mathematical model with reduced mean values of both i-th and j-th input random variables. All other variables were defined in the previous differencing schemes.

Proof of Theorem 2.

Using the simple one-sided backward differencing defined by Equation (3), one can derive mixed partial derivatives as follows:

$$\frac{{\partial }^{2}r}{\partial X_i\partial X_j}\approx \frac{\frac{\partial r\left(X\right)}{\partial X_j}\left({\mu }_{X_i},{\mu }_{X_j}\right)-\frac{\partial r\left(X\right)}{\partial X_j}\left(X_{i\Delta },{\mu }_{X_j}\right)}{{\Delta }_{X_i}}$$

(12)

where $\frac{\partial r\left(X\right)}{\partial X_j}$ is computed for specific coordinates $\left({\mu }_{X_i},{\mu }_{X_j}\right)$ and $\left(X_{i\Delta },{\mu }_{X_j}\right)$ as

$$\frac{\partial r\left(X\right)}{\partial X_j}\left({\mu }_{X_i},{\mu }_{X_j}\right)\approx \frac{R_{X_m}-R_{X_{j\Delta }}}{{\Delta }_{X_j}}$$

(13)

and

$$\frac{\partial r\left(X\right)}{\partial X_j}\left(X_{i\Delta },{\mu }_{X_j}\right)\approx \frac{R_{X_{i\Delta }}-R_{X_{i\Delta }{X}_{j\Delta }}}{{\Delta }_{X_j}}$$

(14)

Therefore, the final derivative scheme for mixed second partial derivatives based on the simple backward differencing adapted for a semi-probabilistic approach is

$$\frac{{\partial }^{2}r\left(X\right)}{\partial X_i\partial X_j}=\frac{R_{X_m}-R_{X_{i\Delta }}-R_{X_{j\Delta }}+R_{X_{i\Delta }{X}_{j\Delta }}}{{\Delta }_{X_i}{\Delta }_{X_j}}$$

(15)

 □

3. Numerical Computation

3.1. Methodology of ECoV by TSE

Since TSE can be constructed in various forms, it is beneficial to create ECoV methodology using TSE, composed of the three levels of an approximation using asymmetric differencing schemes already described in the previous section in combination with linear and quadratic TSE as follows:

(1)

linear TSE with a simple differencing scheme using Equation (3)—$n_{sim}=n+1$,

(2)

linear TSE with an advanced differencing scheme using Equation (9)—$n_{sim}=2n+1$,

(3)

TSE truncated to quadratic terms with a differencing scheme using Equation (9) for the first-order derivatives, Equation (10) for the second-order partial derivatives, and Equation (11) for the mixed derivatives—number of calculation is $n_{sim}=2n+\left(\genfrac{}{}{0pt}{}{n}{2}\right)+1$ in total.

The first level was proposed by Schlune et al. [9] for uncorrelated random variables, and it was used in several practical studies [19,20,21]. However, its behavior for functions of correlated input random variables has not been investigated yet, though it is often necessary to assume correlated random material characteristics in industrial applications. It can be expected that the accuracy of the first level is not sufficient for dependent variables, which will be investigated in numerical examples.

The second level with the advanced differencing scheme still uses only linear terms of the TSE, and thus it is possible to calculate variance by the simple Equation (8), which might be important for easy applications in industry. The accuracy of the second level is significantly improved by additional simulations; however, interaction terms are missing due to a linear truncation of TSE.

The third level of approximation is especially suitable for mathematical models with strong interaction among random variables. However, it is also the most expensive approach, and statistical moments of the model response should be obtained numerically since an analytical calculation is inefficient. Therefore, it can be seen as a simple surrogate model that might be used in combination with Monte Carlo techniques.

Note that the calculations of the original mathematical model from one level are also always used in the following level of approximation. It represents the significant characteristic of the proposed approach, which is beneficial for industrial applications, where it is crucial to decrease the number of calculations as much as possible due to computational demands. Therefore, an analyst can start with the first level of an approximation and eventually increase the number of simulations only if it is necessary. The asymmetric differencing schemes for each level of approximation are depicted in Figure 1 together with iso-lines of bivariate standard Gaussian probability distribution in $\sigma$, $2\sigma$, and $3\sigma$ distance, represented by dotted circles.

3.2. Reference Solution

In industrial applications, only marginal distributions and a correlation matrix are usually known, which does not represent complete information about the joint probability distribution. Therefore, it is necessary to assume a specific copula [22]. A special case of Rosenblatt transformation assuming the Gaussian copula is also known as the Nataf transformation [23], which is usually utilized in reliability applications. The Nataf transformation is composed of three steps:

$$\xi =T_{Nataf}\left(\xi \right)=T_3\circ T_2\circ T_1\left(\xi \right)$$

(16)

The first step represents a transformation from uncorrelated standard Gaussian space $\xi$ to correlated standard normal space $\mathbf{Z}$ using linear transformation.

$$T_1:\xi \mapsto \mathbf{Z}=\mathbf{L}\xi$$

(17)

For this procedure, Cholesky decomposition of the fictive correlation matrix ${\mathbf{R}}_{\mathbf{Z}}$ must be performed:

$${\mathbf{R}}_{\mathbf{Z}}=\mathbf{L}{\mathbf{L}}^T$$

(18)

The following two steps are commonly known as an iso-probabilistic transformation by an inverse cumulative distribution function $F_x^{-1}$ and the standard Gaussian cumulative distribution function $\mathsf{\Phi }$:

$$T_2:\mathbf{Z}\mapsto \mathbf{W}=\mathsf{\Phi }\left(\mathbf{Z}\right)$$

(19)

$$T_3:\mathbf{W}\mapsto \mathbf{X}=F_x^{-1}\left(\mathbf{W}\right)$$

(20)

It is clear that the critical task of the Nataf transformation is to determine ${\mathbf{R}}_{\mathbf{Z}}$. The relationship between the fictive correlation coefficients ${\rho }_{zij}$ and ${\rho }_{ij}$ between i-th and j-th variable is defined by the following integral equation:

$${\rho }_{ij}=\frac{1}{{\sigma }_i{\sigma }_j}{\int \int }_{{\mathbb{R}}^2}\left\{F_i^{-1}\left[\mathsf{\Phi }\left(z_i\right)-{\mu }_i\right]F_j^{-1}\left[\mathsf{\Phi }\left(z_j\right)-{\mu }_j\right]\times {\varphi }_2\left(z_i,z_j,{\rho }_{zij}\right)\right\}d{z}_{i}d{z}_j,$$

(21)

where $\mu$ is the mean value, $\sigma$ is the standard deviation, and ${\varphi }_2$ is the bivariate standard normal probability density function parametrized by fictive correlation coefficients ${\rho }_{zij}$:

$${\varphi }_2(z_i,z_j,{\rho }_{zij})=\frac{1}{2\pi \sqrt{1-{\rho }_{zij}^2}}exp\left(-\frac{z_i^2-2{\rho }_{zij}{z}_i{z}_j+z_j^2}{2(1-{\rho }_{zij}^2)}\right).$$

(22)

Numerical examples are constructed in order to show the behavior of the presented differencing schemes and identify their limitations. For each example, the reference solution is obtained by numerical simulation with $n_{sim}={10}^5$ realizations of a given random vector generated by Latin Hypercube Sampling (LHS) in uncorrelated space $\xi$ and transformed into the correlated space $\mathbf{X}$ by the Nataf transformation. The reference solution by LHS is compared with the results obtained by TSE of increasing complexity using the proposed methodology.

Since this paper is focused on the potential of the presented differencing schemes for industrial applications, the input variables are assumed to be lognormally distributed with coefficient of variation $CoV=0.1–0.2$, which is common for material characteristics. Specifically, all examples work with the following stochastic model of two input variables: vector of mean values $\mu =\left[40,300\right]$ and the corresponding vector of coefficients of variation $CoV=\left[0.1,0.2\right]$. Moreover, Pearson’s correlation coefficients (parameterizing Gaussian copula) are assumed to be positive in the range $\langle 0,0.9\rangle$. The results of the numerical simulations are statistically processed in order to obtain the mean value, variance, and coefficient of variation of the model response.

3.3. Example 1: Simple Linear Model

The very first example represents the entire methodology. It is a simple linear model $R=r\left(\mathbf{X}\right)=X_1+X_2$. The selected realizations generated by Latin Hypercube Sampling, which illustrate the uniform cover of the design domain, together with iso-lines of joint probability density of random vector in uncorrelated and correlated space (Gaussian copula parametrized by the correlation coefficient $\rho =0.8$) are depicted in Figure 2. A reference solution based on a sample with $n_{sim}={10}^5$ is calculated for all examples.

In this case, all the presented differencing schemes led to the exact solution, as can be seen in Figure 3, since a linear approximation fits the original model. For the sake of clarity, the figures in this section show estimation of CoV (top) and variance (bottom) as well, since CoV takes the estimation of mean value into account. The graphs in the right column represent CoV or variance for correlated variables with subtracted uncorrelated values, which represent pure influence of correlation estimated by the presented methods.

3.4. Example 2: Linear Model with Interactions

The second example is focused on the comparison of the first-order and the second-order Taylor series expansion. The first and the second level of approximation use the first-order Taylor expansion; thus, they are not recommended for mathematical models with significant interaction terms, since there are no mixed derivatives in the approximation, and the influence of the interaction is therefore underestimated. The quadratic TSE is the most computationally demanding and the only one reflecting the interaction terms. For the demonstration of this characteristic, the following adaptation of the previous simple mathematical model is assumed:

$$R=r\left(\mathbf{X}\right)=X_1+X_2+5\left(X_1{X}_2\right)$$

(23)

The obtained results are depicted in Figure 4 in the same manner as in the previous example. The estimated mean value for the uncorrelated input random variable was accurate (${\mu }_R=260$). However, using only linear terms of Taylor expansion led to an identical mean value independent of the correlation among input variables. Therefore, the results of CoV are affected by this characteristic, and all methods seem comparable. The accuracy of the used approximations can be clearly seen on the estimation of variance, where the first two levels of an approximation led to identical results, with the error increasing together with the correlation between input random variables. Of course, the obtained results are exact only if the third-level approximation (quadratic Taylor expansion) is used for the estimation of variance, since Hessian of this function is not equal to the zero matrix $\mathbf{0}$. As can be seen, neglecting an interaction among random variables by the first-order Taylor expansion may lead to a significant error in the estimation of statistical moments even for simple linear functions; thus, an analyst should carefully choose the level of approximation in industrial applications considering the nature of the studied physical system.

3.5. Example 3: Approximation of Industrial Example

The third example is motivated by the industrial applications in civil engineering often represented by non-linear finite element models—typically ultimate resistance given by the peak of the load-deflection curve of concrete structural element. The behavior of such a physical system is often monotone with a slightly non-linear progress. A typical function solved by FEM can be found, for example, in [9], and due to the computational demands of FEM, its shape was replicated by the following artificial function:

$$R=r\left(\mathbf{X}\right)=X_1{X}_2-X_1^2-\left(\frac{X_2^2}{30}\right)-\left(X_1-30\right)\left(X_2-200\right)$$

(24)

The exact mean value estimated by LHS was ${\mu }_R=6264$ and by Taylor series ${\mathbb{E}}_{R_T}=6400$, which leads to the difference between the estimation of CoV and variance depicted in Figure 5. However, the estimation of variance and CoV by linear TSE with advanced differencing together with quadratic TSE was accurate. On the other hand, linear Taylor expansion with simple one-sided backward differencing showed a significant error in estimation for all correlation coefficients. The results on the right-hand side of Figure 5 represent the pure influence of correlation, and as can be seen, the slope of the curve estimated by simple linear TSE was significantly different. Thus, this method is not able to correctly identify the role of correlation.

From the previous examples, it is clear that simple linear Taylor expansion as proposed by Schlune et al. is suitable only for functions of uncorrelated variables, which is not a typical industrial problem. However, it is possible to start with simple differencing for uncorrelated problem and add n additional simulations in order to adapt an approximation for correlated variables.

3.6. Example 4: Non-Linear Function

The last example is created in order to show the limitations of all the presented methods with increasing non-linearity of the original mathematical model. The following function has a similar shape as the model in the previous example; however, it is significantly more non-linear:

$$R=r\left(\mathbf{X}\right)=X_1{X}_{2}cos\left(\frac{\pi X_1}{200}\right)cos\left(\frac{\pi X_2}{2000}\right)$$

(25)

The estimated mean value by TSE for the uncorrelated case was ${\mathbb{E}}_{R_T}=8650$, and the exact value estimated by LHS was ${\mu }_R=8468$. Variance and CoV of R estimated by the presented methods are summarized in Figure 6. As can be expected, with higher non-linearity of mathematical models, it was not suitable anymore to use TSE of lower orders as an approximation of the original model. Since computational requirements of higher-order Taylor series expansions are comparable to the commonly known surrogate models, and the estimation of statistical moments is inefficient, one should prefer more advanced surrogate models (e.g., Polynomial Chaos Expansion or Kriging) together with standard statistical methods.

Specifically in this example, the worst results were obtained by the linear TSE with simple differencing, which represents a poor approximation of the original function; thus, the estimation of variance was not satisfied as well. Similarly, a poor accuracy of the estimated influence of correlation can be clearly seen from the different trends of the curves in the column on the right-hand side in Figure 6. However, since there is no significant interaction between the input random variables, the results obtained by linear Taylor series with advanced differencing were almost identical to the more computationally demanding quadratic Taylor series, which might be a crucial advantage in high-dimensional industrial applications solved by FEM.

4. Discussion

The TSE represents a powerful and accurate technique with a strong mathematical background. Unfortunately, it is usually truncated to linear terms in engineering applications, which may generally lead to poor results in the case of non-linear functions and correlated random input variables. Although Schlune et al. proposed the ECoV method based on linear TSE with a simple asymmetric differencing, there are no studies on its limitations and possible generalizations, although TSE is a highly modifiable technique via differencing schemes and a truncation order of an approximation. Therefore, it was necessary to propose different variations of TSE for specific problems and create the novel methodology of three levels of TSE. The proposed methodology was applied on several analytical examples in order to show the limitations of each level. The variations of TSE were proposed with attention to the reduction of computational cost as much as possible, since derivations are computed by finite differencing of FEM in industrial applications. Therefore, each additional level of the methodology works with the information previously obtained from calculations of the original mathematical model; thus, an approximation can be sequentially made more accurate by calculating several additional simulations and combining them with the previous results used in the asymmetric differencing scheme of lower levels of the proposed methodology.

It can be seen from the presented results, that linear TSE fails in the case of significantly non-linear functions (the last example) and functions with important interaction terms (the second example). In such cases, it is necessary to use quadratic TSE (3rd level of proposed methodology) as an approximation. Moreover, the main motivation of this paper is dealing with correlation among random input variables, which has not been investigated yet in the context of ECoV methods. It is clear from the presented examples that the linear TSE with simple differencing (1st level of the proposed methodology) is not suitable for functions of correlated variables. However, once the differencing scheme according to Equation (9) is used in combination with linear TSE (2nd level of the proposed methodology), its accuracy is significantly improved. Thus, if there is not a strong interaction among input random variables, it is not necessary to use quadratic TSE (3rd level of the proposed methodology), which leads to additional computational requirements.

From the point of view of computational costs, it is possible to add higher terms of Taylor series, but it significantly increases the number of derivatives. Therefore, TSE above the second order is inefficient, and advanced surrogate models such as PCE, Krigging, or ANN should be used. On the other hand, a better accuracy of estimation of CoV and variance can be reached by the improved asymmetric differencing scheme as proposed in this paper. Computational requirements are slightly increased from $n+1$, for the traditional scheme according to Equation (3), to $2n+1$ for the proposed scheme according to Equation (9). It is obvious that variance estimation using Equation (9) is significantly improved in comparison to the traditional differencing scheme represented by Equation (3). However, the main advantage of the proposed method is the accuracy of variance estimation in the case of correlated random variables. It is obvious that there is a difference between the curves representing an increment of variance due to correlation (second part of Equation (8)) estimated by both approaches. The difference between both differencing schemes is proportional to a correlation among input random variables; thus, special attention should be given to functions with high correlation among input random variables.

Generally, the proposed methodology proved to be well-suited for typical industrial mathematical models in civil engineering. Moreover, the paper shows the influence of different variants of TSE and level of statistical correlation on estimated CoV, which is a base for semi-probabilistic approaches to determine design value in civil engineering. Such influence can be significant, for more basic random variables certainly amplified, which will be studied on practical examples represented by non-linear finite element models of structures in further research, and the obtained results will be compared to standard normative approaches as PSF and global safety factor method [24] designed specifically for civil engineering.

5. Conclusions

The non-linearity of functions and statistical correlation of input random variables represent crucial aspects in estimating statistical moments of industrial mathematical models. Unfortunately, the accuracy of standard existing methods is not satisfying for such models. Therefore, this paper presents a novel methodology to estimate the coefficient of variation for functions of correlated input random variables. Since mathematical models in civil engineering are often functions of input correlated random variables, it is necessary to develop new and efficient methods based on a semi-probabilistic approach widely accepted for the design and assessment of structures satisfying given safety requirements. Therefore, the methodology of three levels of increasing complexity, accuracy, and computational cost based on Taylor series expansion is proposed and described. The methodology consists of three advanced differencing schemes adapted for civil engineering by step size parameter. The differencing schemes are based on the asymmetric quotient, which is typical for engineering applications, where one is interested in extreme structural behavior leading to failure. The proposed methodology is applied to four analytical examples, and the results are compared to reference solutions obtained by Latin Hypercube Sampling. The analytical examples are constructed in order to show the efficiency and limitations of each differencing scheme: simple linear function, linear function with strong interaction terms, and finally two non-linear functions. From the obtained results, extensively discussed in the previous section, it is clear that it is necessary to choose advanced asymmetric differencing schemes in the cases of correlated input random variables or increase the truncation order of Taylor series expansion. It was shown that its accuracy is significantly higher in comparison to the simple linear TSE (in absolute values but also in a relative trend of influence of correlation). The slight increment of computational demands of the proposed differencing schemes is a significant advantage in comparison to Taylor series of a higher order, where it is necessary to numerically evaluate a large number of additional derivatives. However, it was shown that quadratic TSE is necessary for mathematical models with strong interaction terms.