from Chapter3code import *4 Response Statistics
Utilizing the second-moment information in \(\pmb{\mu}\) and \(\pmb{\Sigma}\), presented in Section 1.11, the objective in this chapter is to calculate the mean and standard deviation of the response in a single structural analysis. The analysis capabilities established in the previous chapters will be employed with the portal frame as a demonstration example. Considering the displacement response, \(u(\mathbf{x})\), where \(\mathbf{x}\) is the vector of random variables, we seek the response statistics \(\mu_u\) and \(\sigma_u\). In doing so, we need to distinguish between linear and nonlinear functions. As explained at the end of Chapter 1, \(u\) is always a nonlinear function of \(\mathbf{x}\) for real structural problems. Again it is highlighted that even in linear structural analysis it is only the load-displacement relationship that is linear.
The fact that \(u(\mathbf{x})\) is a nonlinear function poses a challenge in single-analysis calculation of the response statistics \(\mu_u\) and \(\sigma_u\). Approximate solutions are obtained by first developing a Taylor approximation for \(u(\mathbf{x})\) at the mean, i.e., at \(\mathbf{x}=\pmb{\mu}\). For subsequent derivations it is convenient to switch from vector-matrix notation to index notation; that transition is shown here for the second-order Taylor approximation of the response of interest:
\[ \begin{aligned} u(\mathbf{x}) &\approx u_o + \nabla u_o^{\top} (\mathbf{x}-\pmb{\mu}) + \frac{1}{2} (\mathbf{x}-\pmb{\mu})^{\top} \mathbf{H} (\mathbf{x}-\pmb{\mu}) \\ &= u_o + \frac{\partial u_o}{\partial x_i} (x_i-\mu_i) + \frac{1}{2} (x_i-\mu_i) \frac{\partial^2 u_o}{\partial x_i \partial x_j} (x_j-\mu_j) \end{aligned} \tag{4.1}\]
The notation \(u_o\) highlights that the response and its first- and second-order derivatives are evaluated at the mean. Also notice the summation implied by repeated indices in each term in Equation 4.1. As an example of that rule, the system of equations \(\mathbf{K}\mathbf{u}=\mathbf{F}\) is written \(K_{ij}u_j=F_i\) in index notation, with summation implied over the index \(j\), which appears twice in the left-hand side term. \(j\) is therefore referred to as a dummy index, while \(i\) is the free index that gives the dimension of the result of the multiplication.
Second-moment information for the input variables for the portal frame in Figure 1.2 was given in Listing 3.6. It is now time to utilize that information. The means and standard deviations are readily available from Listing 3.6 and the two functions defined below extract the covariance and correlation matrices:
means, stdvs, distributions, correlation, trackNode, trackDOF, DDMs = linearFrameVariableSpecs()
def getR(means, correlation):
correlation = np.array(correlation)
R = np.identity(len(means))
for i in range(len(correlation)):
rv_i = int(correlation[i, 0])
rv_j = int(correlation[i, 1])
rho = correlation[i, 2]
R[rv_i-1, rv_j-1] = rho
R[rv_j-1, rv_i-1] = rho
return R
def getCovMatrix(means, stdvs, correlation):
if len(correlation) == 0:
covarianceMatrix = np.diag(np.array(stdvs)**2)
else:
D = np.diag(stdvs)
R = getR(means, correlation)
covarianceMatrix = np.dot(np.dot(D, R), D)
return covarianceMatrix
covarianceMatrix = getCovMatrix(means, stdvs, correlation)4.1 Expectation
The statistical moments introduced in Section 1.11 are examples of the more general concept of expectation. In fact, all the moments given as integrals in Table 1.2 and Equation 1.23 can be expressed by the expectation operator \(\mathrm{E}[\,]\):
- Mean: \(\mathrm{E}[x] \equiv \mu\)
- Mean square: \(\mathrm{E}[x^2]\)
- Variance: \(\mathrm{E}[(x-\mu)^2] \equiv \sigma^2 \equiv \mathrm{Var}[x]\)
- Covariance: \(\mathrm{E}[(x_i-\mu_i)(x_j-\mu_j)] \equiv \mathrm{Cov}[x_i, x_j] = \rho_{12}\sigma_i\sigma_j\)
Comparing that list with the integrals in Table 1.2 and Equation 1.23 reveals the definition of the expectation in terms of integration:
\[ \mathrm{E}[(\cdots)] = \int_{-\infty}^{\infty} (\cdots) \cdot f(x) \,dx \tag{4.2}\]
In words, the expectation of a quantity is an integral with the integrand being that quantity multiplied by the PDF of the intervening random variable(s). To that end, the expectations sought in this chapter are the mean,
\[ \mu_u = \mathrm{E}[u(\mathbf{x})] = \int_{-\infty}^{\infty} u(\mathbf{x}) \cdot f(\mathbf{x}) \,dx \tag{4.3}\]
and the variance,
\[ \sigma_u^2 = \mathrm{Var}[u(\mathbf{x})] = \mathrm{E}[(u(\mathbf{x}) - \mu_u)^2] = \int_{-\infty}^{\infty} (u(\mathbf{x}) - \mu_u)^2 \cdot f(\mathbf{x}) \,dx \tag{4.4}\]
where \(f(\mathbf{x})\) is the joint PDF for the random variables. However, in the following, \(f(\mathbf{x})\) is not established and the integration is not carried out. Instead, we employ the Taylor expansion in Equation 4.1, after first deriving five properties of expectation integrals, helpful in the subsequent derivation of \(\mu_u\) and \(\sigma_u^2\):
- The area underneath a PDF is unity:
\[ \mathrm{E}[1] = \int_{-\infty}^{\infty} f(x) \,dx = 1 \tag{4.5}\]
- The expectation of a constant equals that constant:
\[ \mathrm{E}[c] = \int_{-\infty}^{\infty} c \cdot f(x) \,dx = c \tag{4.6}\]
- The expectation of a constant times a random variable is that constant times the mean of the random variable:
\[ \mathrm{E}[c\cdot x] = \int_{-\infty}^{\infty} c\cdot x\cdot f(x) \,dx = c \cdot \mathrm{E}[x] = c \cdot \mu \tag{4.7}\]
- The marginal PDF for a random variable is obtained from a two-variable joint PDF by “integrating out” the other variable:
\[ f(x_i) = \int_{-\infty}^{\infty} x_j \cdot f(x_i, x_j) \,dx_j \tag{4.8}\]
- The covariance between two random variables, \(x_i\) and \(x_j\), is the mean product, \(\mathrm{E}[x_i \cdot x_j]\), minus the product of the means:
\[ \begin{aligned} \mathrm{Cov}[x_i, x_j] &= \mathrm{E}[(x_i-\mu_i)(x_j-\mu_j)] \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (x_i-\mu_i) \cdot (x_j-\mu_j) \cdot f(x_i, x_j) \,dx_i \,dx_j \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} (x_i\cdot x_j - x_i\cdot \mu_j - x_j\cdot \mu_i + \mu_i \cdot \mu_j) \cdot f(x_i, x_j) \,dx_i \,dx_j \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x_i\cdot x_j \cdot f(x_i, x_j) \,dx_i \,dx_j - \mu_j\int_{-\infty}^{\infty} x_i \int_{-\infty}^{\infty} f(x_i, x_j) \,dx_j \,dx_i \\ &- \mu_i\int_{-\infty}^{\infty}x_j \int_{-\infty}^{\infty} f(x_i, x_j) \,dx_i \,dx_j + \mu_i \cdot \mu_j \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x_i, x_j) \,dx_i \,dx_j \\ &= \mathrm{E}[x_i\cdot x_j] - \mu_j\int_{-\infty}^{\infty} x_i \cdot f(x_i) \,dx_i - \mu_i\int_{-\infty}^{\infty}x_j \cdot f(x_j) \,dx_j + \mu_i \cdot \mu_j \cdot 1 \\ &= \mathrm{E}[x_i\cdot x_j] - \mu_i \cdot \mu_j \end{aligned} \tag{4.9}\]
Those five properties are applied to derive the results in the sections below.
4.2 First-order Response Mean
The expectation of the first two terms in the right-hand side of Equation 4.1, i.e., the first-order approximation of the displacement response gives
\[ \begin{aligned} \mu_u &\approx \mathrm{E}[u_o] + \mathrm{E}\left[\frac{\partial u_o}{\partial x_i} (x_i-\mu_i)\right] \\ &= u_o + \frac{\partial u_o}{\partial x_i}(\mathrm{E}[x_i]-\mu_i) \\ &= u_o \end{aligned} \tag{4.10}\]
That means the first-order approximation of the portal frame displacement is the value calulcated in the previous chapter:
input = createLinearFrameInput(*means)
structuralModel = model(input)
u = linearStaticResponse(structuralModel, trackNode, trackDOF)
print(f"First-order approximation of mean response: {u:.4f}")First-order approximation of mean response: 0.01654.3 Second-order Response Mean
Now including all three terms in the right-hand side of Equation 4.1, it is beneficial to first multiply out the parentheses in the last term of that equation before proceeding:
\[ \begin{aligned} u(\mathbf{x}) &\approx u_o + \frac{\partial u_o}{\partial x_i} (x_i-\mu_i) + \frac{1}{2} \frac{\partial^2 u_o}{\partial x_i \partial x_j} (x_i\cdot x_j - x_i\cdot \mu_j - x_j\cdot \mu_i + \mu_i\cdot \mu_j) \end{aligned} \tag{4.11}\]
Taking the expectation yields
\[ \begin{aligned} \mu_u &\approx u_o + \frac{\partial u_o}{\partial x_i} (\mathrm{E}[x_i]-\mu_i) + \frac{1}{2} \frac{\partial^2 u_o}{\partial x_i \partial x_j} \cdot \mathrm{E}[x_i\cdot x_j - x_i\cdot \mu_j - x_j\cdot \mu_i + \mu_i\cdot \mu_j] \\ &= u_o + \frac{1}{2} \frac{\partial^2 u_o}{\partial x_i \partial x_j} \left(\mathrm{E}[x_i\cdot x_j] - \mathrm{E}[x_i]\cdot \mu_j - \mathrm{E}[x_j]\cdot \mu_i + \mu_i\cdot \mu_j \right) \\ &= u_o + \frac{1}{2} \frac{\partial^2 u_o}{\partial x_i \partial x_j} \left(\mathrm{E}[x_i\cdot x_j] - \mu_i\cdot \mu_j \right) \\ &= u_o + \frac{1}{2} \frac{\partial^2 u_o}{\partial x_i \partial x_j} \mathrm{Cov}[x_i, x_j] \\ &= u_o + \frac{1}{2} \cdot \mathbf{H} \odot \pmb{\Sigma} \end{aligned} \tag{4.12}\]
where the Hadamard product, i.e., the sum of component-wise multiplication of the Hessian matrix of second-order response sensitivities and the covariance matrix is identified in the final result. Utilizing the Hessian matrix calculated in the previous chapter, for the portal frame example, the updated mean is
u, void, Hessian = linearStaticSecondOrder(structuralModel, trackNode, trackDOF, DDMs)
meanSO = u + 0.5 * np.sum(Hessian * covarianceMatrix)
print(f"Second-order approximation of mean response: {meanSO:.4f}")Second-order approximation of mean response: 0.01674.4 First-order Response Variance
The inclusion of the last term in Equation 4.1, i.e., the second-order term, is awkward in the calculation of the response variance. Higher statistical moments, rarely available in practice, appear in the result. Therefore, the expression for the response variance is derived using the first two terms in the right-hand side of Equation 4.1:
\[ \begin{aligned} \sigma_u^2 &= \mathrm{Var}[u(\mathbf{x})] = \mathrm{E}[(u(\mathbf{x})-\mu_u)^2] \\ &\approx \mathrm{E}\left[\left(u_o + \frac{\partial u_o}{\partial x_i} (x_i-\mu_i)-\mu_u\right)^2\right] \\ &= \mathrm{E}\left[\left(\frac{\partial u_o}{\partial x_i} (x_i-\mu_i)\right)^2\right] \\ &= \mathrm{E}\left[\left(\frac{\partial u_o}{\partial x_i} x_i-\frac{\partial u_o}{\partial x_i} \mu_i\right)^2\right] \\ &= \mathrm{E}\left[\left(\frac{\partial u_o}{\partial x_i} x_i-\frac{\partial u_o}{\partial x_j} \mu_j\right)\left(\frac{\partial u_o}{\partial x_k} x_k-\frac{\partial u_o}{\partial x_l} \mu_l\right)\right] \\ &= \mathrm{E}\left[\frac{\partial u_o}{\partial x_i} x_i\frac{\partial u_o}{\partial x_k} x_k- \frac{\partial u_o}{\partial x_j} \mu_j\frac{\partial u_o}{\partial x_k} x_k-\frac{\partial u_o}{\partial x_i} x_i\frac{\partial u_o}{\partial x_l} \mu_l+\frac{\partial u_o}{\partial x_j} \mu_j\frac{\partial u_o}{\partial x_l} \mu_l\right] \\ &= \frac{\partial u_o}{\partial x_i} \frac{\partial u_o}{\partial x_k} \mathrm{E}[x_i\cdot x_k]- \frac{\partial u_o}{\partial x_j} \mu_j\frac{\partial u_o}{\partial x_k} \mathrm{E}[x_k]-\frac{\partial u_o}{\partial x_i} \mathrm{E}[x_i]\frac{\partial u_o}{\partial x_l} \mu_l+\frac{\partial u_o}{\partial x_j} \mu_j\frac{\partial u_o}{\partial x_l} \mu_l \\ &= \frac{\partial u_o}{\partial x_i} \frac{\partial u_o}{\partial x_k} \left(\mathrm{E}[x_i\cdot x_k]-\mu_i \mu_k\right) \\ &= \frac{\partial u_o}{\partial x_i} \frac{\partial u_o}{\partial x_k} \cdot \mathrm{Cov}[x_i, x_k] \\ &= \nabla u_o^{\top} \pmb{\Sigma} \nabla u_o \end{aligned} \tag{4.13}\]
In the longest line of that equation, notice that the last three terms are identical, when dummy indices are renamed. The last equality translates the result from index notation into vector-matrix notation, giving the following computational result:
u, dudx = linearStaticFirstOrder(structuralModel, trackNode, trackDOF, DDMs)
stdvFO = np.sqrt(dudx.dot(covarianceMatrix.dot(dudx)))
print(f"First-order response standard deviation: {stdvFO:.4f}")
print(f"({(stdvFO/meanSO*100):.2f}% coefficient of variation)")First-order response standard deviation: 0.0036
(21.52% coefficient of variation)4.5 Sampling
A goal behind this book is to gain insights into the reliability and sensitivity of structures from a single, or just a few, structural analyses. The alternative is brute force sampling, which means repeating the structural analysis many times. Here is a sampling algorithm that produces results that can be compared with the previously calculated mean and standard deviation:
numSamples = int(1e5)
responses = []
H = 6
L = 10
for i in range(numSamples):
Esample = np.random.normal(means[0], stdvs[0])
Asample = np.random.normal(means[1], stdvs[1])
Isample = np.random.normal(means[2], stdvs[2])
qsample = np.random.normal(means[3], stdvs[3])
Fsample = np.random.normal(means[4], stdvs[4])
input = createLinearFrameInput(Esample, Asample, Isample, qsample, Fsample)
structuralModel = model(input)
uSample = linearStaticResponse(structuralModel, trackNode, trackDOF)
responses.append(uSample)
sampling_mean = np.mean(responses)
sampling_std = np.std(responses)
print(f"Sampling mean: {sampling_mean:.4f}")
print(f"Sampling standard deviation: {sampling_std:.4f}")
print(f"({(sampling_std/sampling_mean*100):.2f}% coefficient of variation)")Sampling mean: 0.0167
Sampling standard deviation: 0.0031
(18.29% coefficient of variation)That sampling result confirms that the second-order approximation of the mean response is accurate. Also, it suggests that the first-order approximation of the response standard deviation is a slight overestimation.