![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open"){kind=icon}
![\"Download\"](\"https://img.shields.io/badge/Github-Download-blue.svg\" "\\"Download"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2-Advanced-First-Order-Seccond-Momemt-Methods-(AFSOM))",
"section": "8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)"
}
},
"source": [
"# 8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2-Advanced-First-Order-Seccond-Momemt-Methods-(AFSOM))",
"section": "8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)"
}
},
"source": [
"Created by Qing Lan(qlan@nd.edu), July 3rd, 2020"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2-Advanced-First-Order-Seccond-Momemt-Methods-(AFSOM))",
"section": "8.2 Advanced First-Order Seccond-Momemt Methods (AFSOM)"
}
},
"source": [
"The following text, example, and code were adapted from:\n",
"\n",
"McClarren, Ryan G (2018). *Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers*, *Chapter 8: Reliability Methods for Estimating the Probability of Failure*, Springer, https://doi.org/10.1007/978-3-319-99525-0_8"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"source": [
"## 8.2.1 Reproduce Figure 8.6: "
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"source": [
"$Q(\\mathbf{X}) = 2x_1^3+10x_1 x_2 + x_1 + 3x_2^3+x_2$ \n",
"\n",
"$x_1 \\backsim \\mathcal{N} (0.1, 2^2)$ \n",
"\n",
"$x_2 \\backsim \\mathcal{N} (-0.05, 3^2)$ \n",
"\n",
"Covariance matrix $\\Sigma = \\begin{pmatrix}\n",
"4 & 3.9\\\\\n",
"3.9 & 9\n",
"\\end{pmatrix}$\n",
"\n",
"Correlation matrix $ \\mathbf{R} = \\begin{pmatrix}\n",
"1 & 0.65\\\\\n",
"0.65 & 1\n",
"\\end{pmatrix}$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.stats import norm"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"outputs": [],
"source": [
"# Define the function for QoI\n",
"def QoI(x):\n",
" x1 = x[0].item()\n",
" x2 = x[1].item()\n",
" return 2*x1**3 + 10*x1*x2 + x1 + 3*x2**3 + x2 # a scalar\n",
"\n",
"# Derivatives of QoI with respect to X\n",
"def DQ_DX(x):\n",
" x1 = x[0].item()\n",
" x2 = x[1].item()\n",
" dQ_dx1 = 6*x1**2 + 10*x2 + 1\n",
" dQ_dx2 = 10*x1 + 9*x2**2 + 1\n",
" return np.array([[dQ_dx1, dQ_dx2]]) # a 2-D column array\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"source": [
"Standardize the coordinate system:\n",
"\n",
"$Y_i = \\frac{x_i-\\mu_i'}{\\sigma_i'}$\n",
"\n",
"so that\n",
"\n",
"$Y_i \\backsim \\mathcal{N} (0, 1)$\n",
"\n",
"with known correlation matrix $\\mathbf{R}$ (unchanged)\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"outputs": [],
"source": [
"# Derivatives of QoI to Y using chain rule\n",
"def DQ_DY(x):\n",
" x1 = x[0].item()\n",
" x2 = x[1].item()\n",
" dQ_dy1 = 12*x1**2 + 20*x2 + 2\n",
" dQ_dy2 = 30*x1 + 27*x2**2 + 3\n",
" return np.array([[dQ_dy1, dQ_dy2]]) # a 2-D column array"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"source": [
"Failure surface:\n",
"\n",
"$Z = Q_{fail}-Q{(\\mathbf{X})} = 0$ \n",
"\n",
"We seek for finding the nearest point on the failure surface to the nominal value, which means minimizing\n",
"\n",
"$\\beta = \\underset{Z(\\mathbf{X})=0}{\\text{min}} \\sqrt{\\mathbf{Y^T(X)R^{-1}Y(X)}}$\n",
"\n",
"With Lagrange multiplier method in Eq.8.11, we get\n",
"\n",
"$\\beta = \\frac{Q_{fail}-Q(\\mathbf{X})+\\nabla_\\mathbf{Y}QY} {\\sqrt{\\nabla_\\mathbf{Y}Q \\mathbf{R} \\nabla_\\mathbf{Y}^{T}Q}}$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"nbpages": {
"level": 2,
"link": "[8.2.1 Reproduce Figure 8.6: ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/08.02-Advanced-First-Order-Second-Moment-Method.html#8.2.1-Reproduce-Figure-8.6:)",
"section": "8.2.1 Reproduce Figure 8.6: "
}
},
"outputs": [],
"source": [
"def AFOSM(x0, pdf_x, cdf_x, R, Q_fail, delta, eps, max_steps):\n",
" \n",
" '''\n",
" x0: initial guess\n",
" pdf_x: list of margin PDF functions of variables in x\n",
" cdf_x: list of margin CDF functions of variables in x\n",
" R: correlation matrix of x\n",
" Q_fail: failure threshold\n",
" delta: tolerance of convergence of beta\n",
" eps: tolerance of convergence of abs(Q-Q_fail)\n",
" \n",
" '''\n",
" \n",
" \n",
" \n",
" def solv_y (xl, pdf_x, cdf_x):\n",
" # a function that finds standardized variable y and the corresponding mu_y and sigma_y\n",
" # based on the fact that CDF and PDF at ccorresponding points should be the same as the originals\n",
" \n",
" '''\n",
" inputs: x, mu_x, sigma_x\n",
" outputs: yi=(xi-mu_xi), mu_y, sigma_y\n",
" '''\n",
" \n",
" # get dimension of inputs\n",
" dim = len(x0)\n",
" \n",
" # initialize new variable Y, mu_y and sigma_y\n",
" yl = x0*0\n",
" sigma_y = x0*0\n",
" mu_y = x0*0\n",
" \n",
" # solve for mu_y and sigma_y\n",
" for i in range(dim):\n",
" xi = xl[i].item()\n",
" \n",
" # get the distribution functions from the list\n",
" fun_pdf = pdf_x[i]\n",
" fun_cdf = cdf_x[i]\n",
" # cumulative distribution function value at x_i\n",
" cdf_xi = fun_cdf(xi)\n",
" # probability distribution function value at x_i\n",
" pdf_xi = fun_pdf(xi)\n",
" \n",
" # solution for mu and sigma of yi, using the formula Eq.8.8 and 8.9\n",
" sigma_yi = norm.pdf(norm.ppf(cdf_xi))/pdf_xi\n",
" mu_yi = xi-norm.ppf(cdf_xi)*sigma_yi\n",
"\n",
" # save the solution into arrays\n",
" sigma_y[i] = sigma_yi \n",
" mu_y[i] = mu_yi\n",
" # standardized coordinate\n",
" yl[i] = (xi-mu_yi)/sigma_yi\n",
" \n",
" return yl, mu_y, sigma_y\n",
" \n",
"\n",
" \n",
" # start up with the first two steps\n",
" \n",
" # collect steps\n",
" beta = []\n",
" xl_lst = []\n",
" # initial step\n",
" l=0\n",
" xl = x0\n",
" xl_lst.append(x0)\n",
" # solve for yl from the current xl\n",
" yl, mu_y, sigma_y = solv_y (xl, pdf_x, cdf_x)\n",
" beta1 = np.sqrt((yl.T@np.linalg.inv(R)@yl).item())\n",
" # Lagrange mutiplier\n",
" lamda = (Q_fail - QoI(xl) + (DQ_DY(xl)@yl).item())/(DQ_DY(xl)@R@DQ_DY(xl).T).item()\n",
" beta.append(beta1)\n",
" \n",
" # second step\n",
" l = l+1\n",
" # update yl, beta2\n",
" yl = lamda*R@(DQ_DY(xl).T)\n",
" beta2 = np.sqrt((yl.T@np.linalg.inv(R)@yl).item())\n",
" beta.append(beta2)\n",
" # update xl\n",
" xl = sigma_y*yl + mu_y\n",
" xl_lst.append(xl)\n",
" \n",
" # iterate until convergence\n",
" while l"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}