Merge pull request #225 from SURGroup/feature/inverse_form

Feature/inverse form

docs/code/reliability/inverse_form/README.rst

@@ -0,0 +1,2 @@
+Inverse FORM Examples
+^^^^^^^^^^^^^^^^^^^^^^^^^

docs/code/reliability/inverse_form/inverse_form_cantilever.py

@@ -0,0 +1,73 @@
+"""
+Inverse FORM - Cantilever Beam
+-----------------------------------
+The following example is Example 7.2 from Chapter 7 of :cite:`FORM_XDu`.
+
+A cantilever beam is considered to fail if the displacement at the tip exceeds the threshold :math:`D_0`.
+The performance function :math:`G(\\textbf{U})` of this problem is given by
+
+.. math:: G = D_0 - \\frac{4L^3}{Ewt} \\sqrt{ \\left(\\frac{P_x}{w^2}\\right)^2 + \\left(\\frac{P_y}{t^2}\\right)^2}
+
+Where the external forces are modeled as random variables :math:`P_x \sim N(500, 100)` and :math:`P_y \sim N(1000,100)`.
+The constants in the problem are length (:math:`L=100`), elastic modulus (:math:`E=30\\times 10^6`), cross section width
+(:math:`w=2`), cross section height (:math:`t=4`), and :math:`D_0=3`.
+
+"""
+# %% md
+#
+# First, we import the necessary modules.
+
+# %%
+
+import numpy as np
+from scipy import stats
+from UQpy.distributions import Normal
+from UQpy.reliability.taylor_series import InverseFORM
+from UQpy.run_model.RunModel import RunModel
+from UQpy.run_model.model_execution.PythonModel import PythonModel
+
+# %% md
+#
+# Next, we initialize the :code:`RunModel` object.
+# The file defining the performance function file can be found on the UQpy GitHub.
+# It contains a function :code:`cantilever_beam` to compute the performance function :math:`G(\textbf{U})`.
+
+# %%
+
+model = PythonModel(model_script='local_pfn.py', model_object_name="cantilever_beam")
+runmodel_object = RunModel(model=model)
+
+# %% md
+#
+# Next, we define the external forces in the :math:`x` and :math:`y` direction as distributions that will be passed into
+# :code:`FORM`. Along with the distributions, :code:`FORM` takes in the previously defined :code:`runmodel_object`,
+# the specified probability of failure, and the tolerances. These tolerances are smaller than the defaults to ensure
+# convergence with the level of accuracy given in the problem.
+
+# %%
+
+p_fail = 0.04054
+distributions = [Normal(500, 100), Normal(1_000, 100)]
+inverse_form = InverseFORM(distributions=distributions,
+                           runmodel_object=runmodel_object,
+                           p_fail=p_fail,
+                           tolerance_u=1e-5,
+                           tolerance_gradient=1e-5)
+
+# %% md
+#
+# With everything defined we are ready to run the inverse first-order reliability method and print the results.
+# The solution to this problem given by Du is :math:`\textbf{U}^*=(1.7367, 0.16376)` with a reliability index of
+# :math:`\beta_{HL}=||\textbf{U}^*||=1.7444` and probability of failure of
+# :math:`p_{fail} = \Phi(-\beta_{HL})=0.04054`. We expect this problem to converge in 4 iterations.
+# We confirm our design point matches this length, and therefore has a probability of failure specified by our input.
+
+# %%
+
+inverse_form.run()
+beta = np.linalg.norm(inverse_form.design_point_u)
+print('Design point in standard normal space (u^*):', inverse_form.design_point_u[0])
+print('Design point in original space:', inverse_form.design_point_x[0])
+print('Hasofer-Lind reliability index:', beta)
+print('Probability of failure at design point:', stats.norm.cdf(-beta))
+print('Number of iterations:', inverse_form.iteration_record[0])

docs/code/reliability/inverse_form/local_pfn.py

@@ -0,0 +1,24 @@
+"""
+
+Auxiliary file
+==============================================
+
+"""
+import numpy as np
+
+
+def cantilever_beam(samples=None):
+    """Performance function from Chapter 7 Example 7.2 from Du 2005"""
+    elastic_modulus = 30e6
+    length = 100
+    width = 2
+    height = 4
+    d_0 = 3
+
+    g = np.zeros(samples.shape[0])
+    for i in range(samples.shape[0]):
+        x = (samples[i, 0] / width**2) ** 2
+        y = (samples[i, 1] / height**2) ** 2
+        d = ((4 * length**3) / (elastic_modulus * width * height)) * np.sqrt(x + y)
+        g[i] = d_0 - d
+    return g

docs/source/conf.py

@@ -106,6 +106,7 @@
         "../code/reliability/form",
         "../code/reliability/sorm",
         "../code/reliability/subset_simulation",
+        "../code/reliability/inverse_form",
         "../code/surrogates/srom",
         "../code/surrogates/gpr",
         "../code/surrogates/pce",
@@ -148,6 +149,7 @@
         "auto_examples/reliability/form",
         "auto_examples/reliability/sorm",
         "auto_examples/reliability/subset_simulation",
+        "auto_examples/reliability/inverse_form",
         "auto_examples/surrogates/srom",
         "auto_examples/surrogates/gpr",
         "auto_examples/surrogates/pce",

docs/source/reliability/index.rst

@@ -2,10 +2,10 @@ Reliability
 ===========
 
 Reliability of a system refers to the assessment of its probability of failure (i.e the system no longer satisfies some
-performance measure), given the model uncertainty in the structural, environmental and load parameters. Given a vector
+performance measures), given the model uncertainty in the structural, environmental and load parameters. Given a vector
 of random variables :math:`\textbf{X}=[X_1, X_2, \ldots, X_n]^T \in \mathcal{D}_\textbf{X}\subset \mathbb{R}^n`, where
-:math:`\mathcal{D}_\textbf{X}` is the domain of interest and :math:`f_{\textbf{X}}(\textbf{x})` is its joint probability density
-function, then the probability that the system will fail is defined as
+:math:`\mathcal{D}` is the domain of interest and :math:`f_{\textbf{X}}(\textbf{x})` is its joint probability density
+function then, the probability that the system will fail is defined as
 
 
 .. math:: P_f =\mathbb{P}(g(\textbf{X}) \leq 0) = \int_{\mathcal{D}_f} f_{\textbf{X}}(\textbf{x})d\textbf{x} = \int_{\{\textbf{X}:g(\textbf{X})\leq 0 \}} f_{\textbf{X}}(\textbf{x})d\textbf{x}
@@ -23,8 +23,8 @@ transformation and can support reliability analysis for problems with arbitraril
 This module contains functionality for all reliability methods supported in :py:mod:`UQpy`.
 The module currently contains the following classes:
 
-- :class:`.TaylorSeries`: Class to perform reliability analysis using First Order reliability Method (:class:`FORM`) and Second Order
-  Reliability Method (:class:`SORM`).
+- :class:`.TaylorSeries`: Class to perform reliability analysis using First Order Reliability Method (:class:`FORM`),
+  Inverse First Order Reliability Method (:class:`InverseFORM`) and Second Order Reliability Method (:class:`SORM`).
 - :class:`.SubsetSimulation`: Class to perform reliability analysis using subset simulation.
 
 

docs/source/reliability/inverse_form.rst

@@ -0,0 +1,77 @@
+Inverse FORM
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+In FORM :cite:`FORM_XDu`, the performance function is linearized according to
+
+.. math:: G(\textbf{U})  \approx  G(\textbf{U}^\star) + \nabla G(\textbf{U}^\star)(\textbf{U}-\textbf{U}^\star)^\intercal
+
+where :math:`\textbf{U}^\star` is the expansion point, :math:`G(\textbf{U})` is the performance function evaluated in
+the standard normal space and :math:`\nabla G(\textbf{U}^\star)` is the gradient of :math:`G(\textbf{U})` evaluated at
+:math:`\textbf{U}^\star`. The probability failure is approximated as
+
+.. math:: p_{fail} = \Phi(-\beta_{HL})
+
+where :math:`\Phi(\cdot)` is the standard normal cumulative distribution function and :math:`\beta_{HL}=||\textbf{U}^*||`
+is the norm of the design point known as the Hasofer-Lind reliability index.
+
+The goal of the inverse FORM algorithm is to find a design point :math:`\textbf{U}^\star` that minimizes the performance
+function :math:`G(\textbf{U})` and lies on the hypersphere defined by :math:`||\textbf{U}^*|| = \beta_{HL}`, or
+equivalently :math:`||\textbf{U}^*|| = -\Phi^{-1}(p_{fail})`. The default convergence criteria used for this algorithm
+are:
+
+.. math:: \text{tolerance}_{\textbf{U}}:\ ||\textbf{U}_i - \textbf{U}_{i-1}||_2 \leq 10^{-3}
+.. math:: \text{tolerance}_{\nabla G(\textbf{U})}:\ ||\nabla G(\textbf{U}_i)- \nabla G(\textbf{U}_{i-1})||_2 \leq 10^{-3}
+
+
+**Problem Statement**
+
+Compute :math:`u^* = \text{argmin}\ G(\textbf{U})` such that :math:`||\textbf{U}||=\beta`.
+
+The feasibility criteria :math:`||\textbf{U}||=\beta` may be equivalently defined as
+:math:`\beta = -\Phi^{-1}(p_{fail})`, where :math:`\Phi^{-1}(\cdot)` is the inverse standard normal CDF.
+
+**Algorithm**
+
+This method implements a gradient descent algorithm to solve the optimization problem within the tolerances specified by
+:code:`tolerance_u` and :code:`tolerance_gradient`.
+
+0. Define :math:`u_0` and :math:`\beta` directly or :math:`\beta=-\Phi^{-1}(p_\text{fail})`
+1. Set :math:`u=u_0` and :math:`\text{converged}=False`
+2. While not :math:`\text{converged}`:
+    a. :math:`\alpha = \frac{\nabla G(u)}{|| \nabla G(u) ||}`
+    b. :math:`u_{new} = - \alpha \beta`
+    c. If :math:`||u_{new} - u || \leq \text{tolerance}_u` and/or :math:`|| \nabla G(u_{new}) - \nabla G(u) || \leq \text{tolerance}_{\nabla G(u)}`;
+        set :math:`\text{converged}=True`
+       else;
+        :math:`u = u_{new}`
+
+The :class:`.InverseFORM` class is imported using the following command:
+
+>>> from UQpy.reliability.taylor_series import InverseFORM
+
+Methods
+-------
+
+.. autoclass:: UQpy.reliability.taylor_series.InverseFORM
+    :members: run
+
+Attributes
+----------
+
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.alpha_record
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.beta_record
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.design_point_u
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.design_point_x
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.error_record
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.iteration_record
+.. autoattribute:: UQpy.reliability.taylor_series.InverseFORM.failure_probability_record
+
+
+Examples
+--------
+
+.. toctree::
+
+   InverseFORM Examples <../auto_examples/reliability/inverse_form/index>

docs/source/reliability/taylor_series.rst

@@ -1,27 +1,31 @@
 Taylor Series
 -------------
 
-:class:`.TaylorSeries` is a class that calculates the reliability  of a model using the First Order Reliability Method (FORM)
-or the Second Order Reliability Method (SORM) based on the first-order and second-order Taylor series expansion
-approximation of the performance function, respectively (:cite:`TaylorSeries1`, :cite:`TaylorSeries2`).
+:class:`.TaylorSeries` is a class that calculates the reliability  of a model using the First Order Reliability Method
+(FORM), Inverse First Order Reliability Method (InverseFORM) or Second Order Reliability Method (SORM) based on the
+first-order and second-order Taylor series expansion approximation of the performance function, respectively
+(:cite:`TaylorSeries1`, :cite:`FORM_XDu`, :cite:`TaylorSeries2`).
 
 .. image:: ../_static/Reliability_FORM.png
    :scale: 40 %
    :alt:  Graphical representation of the FORM.
    :align: center
 
-The :class:`.TaylorSeries` class is the parent class of the :class:`.FORM` and :class:`.SORM` classes that perform the FORM and SORM,
-respectively. These classes can be imported in a python script using the following command:
+The :class:`.TaylorSeries` class is the parent class of the :class:`.FORM`, :class:`InverseFORM`, and :class:`.SORM`
+classes that perform the FORM, InverseFORM, and SORM, respectively.
+These classes can be imported in a python script using the following command:
 
->>> from UQpy.reliability.taylor_series import FORM, SORM
+>>> from UQpy.reliability.taylor_series import FORM, InverseFORM, SORM
 
 
 .. toctree::
    :maxdepth: 1
 
     FORM <form>
+    InverseFORM <inverse_form>
     SORM <sorm>
 
 
 
 
+

src/UQpy/reliability/taylor_series/FORM.py

@@ -133,6 +133,8 @@ def __init__(
         self.x_record: list = []
         """Record of all iteration points in the parameter space **X**."""
 
+        if (seed_x is not None) and (seed_u is not None):
+            raise ValueError('UQpy: Only one input (seed_x or seed_u) may be provided')
         if self.seed_u is not None:
             self.run(seed_u=self.seed_u)
         elif self.seed_x is not None: