{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Probability estimation with importance sampling simulation on cantilever beam example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example we estimate a failure probability with the importance sampling simulation algorithm provided by the `ImportanceSamplingExperiment` class. \n", "\n", "The main steps of this method are:\n", "\n", "* run a FORM analysis,\n", "* create an importance distribution based on the results of the FORM results,\n", "* run a sampling-based probability estimate algorithm. \n", "\n", "\n", "## Description of the problem\n", "\n", "Let us consider the analytical example of a cantilever beam, with Young's modulus E, length L and section modulus I.\n", "\n", "One end of the cantilever beam is built in a wall and we apply a concentrated bending load F at the other end of the beam, resulting in a deviation:\n", "\n", "$$\n", "d = \\frac{FL^3}{3EI}\n", "$$\n", "\n", "Failure occurs when the beam deviation is too large:\n", "\n", "$$d \\ge 30 (cm) $$\n", "\n", "Four independent random variables are considered:\n", "\n", " - E: Young's modulus [Pa]\n", " - F: load [N]\n", " - L: length [m]\n", " - I: section [m^4]\n", "\n", "Stochastic model (simplified model, no units):\n", "\n", " - E ~ Beta(0.93, 2.27, 2.8e7, 4.8e7)\n", " - F ~ LogNormal(30000, 9000, 15000)\n", " - L ~ Uniform(250, 260)\n", " - I ~ Beta(2.5, 1.5, 3.1e2, 4.5e2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define the cantilever beam model" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from __future__ import print_function\n", "import openturns as ot" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Create the marginal distributions of the parameters\n", "dist_E = ot.Beta(0.93, 2.27, 2.8e7, 4.8e7)\n", "dist_F = ot.LogNormalMuSigma(30000, 9000, 15000).getDistribution()\n", "dist_L = ot.Uniform(250, 260)\n", "dist_I = ot.Beta(2.5, 1.5, 3.1e2, 4.5e2)\n", "marginals = [dist_E, dist_F, dist_L, dist_I]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Create the Copula\n", "RS = ot.CorrelationMatrix(4)\n", "RS[2, 3] = -0.2\n", "# Evaluate the correlation matrix of the Normal copula from RS\n", "R = ot.NormalCopula.GetCorrelationFromSpearmanCorrelation(RS)\n", "# Create the Normal copula parametrized by R\n", "copula = ot.NormalCopula(R) " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Create the joint probability distribution\n", "distribution = ot.ComposedDistribution(marginals, copula)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# create the model\n", "model = ot.SymbolicFunction(['E', 'F', 'L', 'I'], ['F*L^3/(3*E*I)'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define the event" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create the event whose probability we want to estimate." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "vect = ot.RandomVector(distribution)\n", "G = ot.CompositeRandomVector(model, vect)\n", "event = ot.ThresholdEvent(G, ot.Greater(), 30.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run a FORM analysis" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Define a solver\n", "optimAlgo = ot.Cobyla()\n", "optimAlgo.setMaximumEvaluationNumber(1000)\n", "optimAlgo.setMaximumAbsoluteError(1.0e-10)\n", "optimAlgo.setMaximumRelativeError(1.0e-10)\n", "optimAlgo.setMaximumResidualError(1.0e-10)\n", "optimAlgo.setMaximumConstraintError(1.0e-10)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Run FORM\n", "algo = ot.FORM(optimAlgo, event, distribution.getMean())\n", "algo.run()\n", "result = algo.getResult()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configure an importance sampling algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `getStandardSpaceDesignPoint` method returns the design point in the U-space." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
[-0.602386,2.31056,0.355794,-0.533677]
" ], "text/plain": [ "class=Point name=Standard Space Design Point dimension=4 values=[-0.602386,2.31056,0.355794,-0.533677]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()\n", "standardSpaceDesignPoint" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The key point is to define the importance distribution in the U-space. To define it, we use a multivariate standard Gaussian centered around the design point in the U-space." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dimension = distribution.getDimension()\n", "dimension" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Normal(mu = [-0.602386,2.31056,0.355794,-0.533677], sigma = [1,1,1,1], R = [[ 1 0 0 0 ]
\n",
" [ 0 1 0 0 ]
\n",
" [ 0 0 1 0 ]
\n",
" [ 0 0 0 1 ]])