Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/inverse form #225

Merged
merged 18 commits into from
Sep 1, 2023
Merged
Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
9cd9d75
initial commit of the InverseFORM class
connor-krill Aug 15, 2023
c3e5459
tests for accuracy and behavior under various inputs for InverseFORM …
connor-krill Aug 15, 2023
48b6204
documentation for InverseFORM as a child of the TaylorSeries class
connor-krill Aug 15, 2023
55db785
scripts to config and run examples for InverseFORM
connor-krill Aug 15, 2023
bd5c57d
added InverseFORM to paths so the examples would show up in the ReadT…
connor-krill Aug 15, 2023
935d45d
improved docstrings and markdown for cleaniness and consistency in re…
connor-krill Aug 15, 2023
b74b418
corrected vector notation in math statement
connor-krill Aug 16, 2023
770256b
edited markdown for clarity and notational consistency
connor-krill Aug 16, 2023
fa02302
editted docstrings for clarity and notational consistency
connor-krill Aug 16, 2023
adafcdf
added references to the new child class InverseFORM
connor-krill Aug 16, 2023
5a6a1ce
added state_function attribute, correct attribute behavior when max n…
connor-krill Aug 17, 2023
d992388
fixed minor bug with state_function attribute
connor-krill Aug 17, 2023
738e25b
changed title font to match existing documentation
connor-krill Aug 18, 2023
6236297
added logic to prevent seed_x and see_u from being specified at the s…
connor-krill Sep 1, 2023
9704aa0
added logic to run InverseFORM.run() if seed is provided in init
connor-krill Sep 1, 2023
1aa56e0
fixes dummy syntax
connor-krill Sep 1, 2023
d8c3b73
deleted unneeded attribute
connor-krill Sep 1, 2023
bc50b1c
Merge remote-tracking branch 'origin/Development' into feature/invers…
dimtsap Sep 1, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/code/reliability/inverse_form/README.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Inverse FORM Examples
^^^^^^^^^^^^^^^^^^^^^^^^^
73 changes: 73 additions & 0 deletions docs/code/reliability/inverse_form/inverse_form_cantilever.py
Original file line number Diff line number Diff line change
@@ -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])
24 changes: 24 additions & 0 deletions docs/code/reliability/inverse_form/local_pfn.py
Original file line number Diff line number Diff line change
@@ -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
2 changes: 2 additions & 0 deletions docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down
6 changes: 3 additions & 3 deletions docs/source/reliability/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ Reliability

Reliability of a system refers to the assessment of its probability of failure (i.e the system no longer satisfies some
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\} \in \mathcal{D}_\textbf{X}\subset \mathbb{R}^n`, where
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}` 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

Expand All @@ -22,8 +22,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.


Expand Down
77 changes: 77 additions & 0 deletions docs/source/reliability/inverse_form.rst
Original file line number Diff line number Diff line change
@@ -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>
16 changes: 10 additions & 6 deletions docs/source/reliability/taylor_series.rst
Original file line number Diff line number Diff line change
@@ -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>





2 changes: 2 additions & 0 deletions src/UQpy/reliability/taylor_series/FORM.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Loading