diff --git a/docs/code/reliability/inverse_form/README.rst b/docs/code/reliability/inverse_form/README.rst new file mode 100644 index 00000000..44cdc4d1 --- /dev/null +++ b/docs/code/reliability/inverse_form/README.rst @@ -0,0 +1,2 @@ +Inverse FORM Examples +^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/code/reliability/inverse_form/inverse_form_cantilever.py b/docs/code/reliability/inverse_form/inverse_form_cantilever.py new file mode 100644 index 00000000..54f355e7 --- /dev/null +++ b/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]) diff --git a/docs/code/reliability/inverse_form/local_pfn.py b/docs/code/reliability/inverse_form/local_pfn.py new file mode 100644 index 00000000..4f40a8a7 --- /dev/null +++ b/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 diff --git a/docs/source/conf.py b/docs/source/conf.py index b4891bf8..2f4e7635 100644 --- a/docs/source/conf.py +++ b/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", diff --git a/docs/source/reliability/index.rst b/docs/source/reliability/index.rst index 51362cfb..1a6666cd 100644 --- a/docs/source/reliability/index.rst +++ b/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. diff --git a/docs/source/reliability/inverse_form.rst b/docs/source/reliability/inverse_form.rst new file mode 100644 index 00000000..502b3c9c --- /dev/null +++ b/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> \ No newline at end of file diff --git a/docs/source/reliability/taylor_series.rst b/docs/source/reliability/taylor_series.rst index 29ff80e0..78ec77a6 100644 --- a/docs/source/reliability/taylor_series.rst +++ b/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
+ InverseFORM SORM + diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index acd94e10..8c167e52 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/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: diff --git a/src/UQpy/reliability/taylor_series/InverseFORM.py b/src/UQpy/reliability/taylor_series/InverseFORM.py new file mode 100644 index 00000000..1e4f4fa9 --- /dev/null +++ b/src/UQpy/reliability/taylor_series/InverseFORM.py @@ -0,0 +1,215 @@ +import logging +import warnings +import numpy as np +import scipy.stats as stats +from typing import Union +from beartype import beartype +from UQpy.distributions import * +from UQpy.transformations import * +from UQpy.run_model.RunModel import RunModel +from UQpy.utilities.ValidationTypes import PositiveInteger +from UQpy.reliability.taylor_series.baseclass.TaylorSeries import TaylorSeries + +warnings.filterwarnings('ignore') + + +class InverseFORM(TaylorSeries): + + @beartype + def __init__( + self, + distributions: Union[None, Distribution, list[Distribution]], + runmodel_object: RunModel, + p_fail: Union[None, float] = 0.05, + beta: Union[None, float] = None, + seed_x: Union[list, np.ndarray] = None, + seed_u: Union[list, np.ndarray] = None, + df_step: Union[int, float] = 0.01, + corr_x: Union[list, np.ndarray] = None, + corr_z: Union[list, np.ndarray] = None, + max_iterations: PositiveInteger = 100, + tolerance_u: Union[float, int, None] = 1e-3, + tolerance_gradient: Union[float, int, None] = 1e-3, + ): + """Class to perform the Inverse First Order Reliability Method. + + Each time :meth:`run` is called the results are appended to the class attributes. + By default, :code:`tolerance_u` and :code:`tolerance_gradient` must be satisfied for convergence. + Specifying a tolerance as :code:`None` will ignore that criteria, but both cannot be :code:`None`. + This is a child class of the :class:`TaylorSeries` class. + + :param distributions: Marginal probability distributions of each random variable. Must be of + type :class:`.DistributionContinuous1D` or :class:`.JointIndependent`. + :param runmodel_object: The computational model. Must be of type :class:`RunModel`. + :param p_fail: Probability of failure defining the feasibility criteria as :math:`||u||=-\\Phi^{-1}(p_{fail})`. + Default: :math:`0.05`. Set to :code:`None` to use :code:`beta` as the feasibility criteria. + Only one of :code:`p_fail` or :code:`beta` may be provided. + :param beta: Hasofer-Lind reliability index defining the feasibility criteria as :math:`||u|| = \\beta`. + Default: :code:`None`. + Only one of :code:`p_fail` or :code:`beta` may be provided. + :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If either :code:`seed_u` or :code:`seed_x` is provided, then the :py:meth:`run` method will be executed automatically. + :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If either :code:`seed_u` or :code:`seed_x` is provided, then the :py:meth:`run` method will be executed automatically. + :param df_step: Finite difference step in standard normal space. Default: :math:`0.01` + :param corr_x: Correlation matrix :math:`\mathbf{C_X}` of the random vector :math:`\mathbf{X}` + :param corr_z: Correlation matrix :math:`\mathbf{C_Z}` of the random vector :math:`\mathbf{Z}` + Default: :code:`corr_z` is the identity matrix. + :param max_iterations: Maximum number of iterations for the `HLRF` algorithm. Default: :math:`100` + :param tolerance_u: Convergence threshold for criterion :math:`||\mathbf{U}_i - \mathbf{U}_{i-1}||_2 \leq` + :code:`tolerance_u` of the `HLRF` algorithm. + Default: :math:`1.0e-3` + :param tolerance_gradient: Convergence threshold for criterion + :math:`||\\nabla G(\mathbf{U}_i)- \\nabla G(\mathbf{U}_{i-1})||_2 \leq` :code:`tolerance_gradient` + of the `HLRF` algorithm. + Default: :math:`1.0e-3` + """ + self.distributions = distributions + self.runmodel_object = runmodel_object + if (p_fail is not None) and (beta is not None): + raise ValueError('UQpy: Exactly one input (p_fail or beta) must be provided') + elif (p_fail is None) and (beta is None): + raise ValueError('UQpy: Exactly one input (p_fail or beta) must be provided') + elif p_fail is not None: + self.p_fail = p_fail + self.beta = -stats.norm.ppf(self.p_fail) + elif beta is not None: + self.p_fail = stats.norm.cdf(-1*beta) + self.beta = beta + self.seed_x = seed_x + self.seed_u = seed_u + self.df_step = df_step + self.corr_x = corr_x + self.corr_z = corr_z + self.max_iterations = max_iterations + self.tolerance_u = tolerance_u + self.tolerance_gradient = tolerance_gradient + if (self.tolerance_u is None) and (self.tolerance_gradient is None): + raise ValueError('UQpy: At least one tolerance (tolerance_u or tolerance_gradient) must be provided') + + self.logger = logging.getLogger(__name__) + self.nataf_object = Nataf(distributions=distributions, corr_z=corr_z, corr_x=corr_x) + + # Determine the number of dimensions as the number of random variables + if isinstance(distributions, DistributionContinuous1D): + self.dimension = 1 + elif isinstance(distributions, JointIndependent): + self.dimension = len(distributions.marginals) + elif isinstance(distributions, list): + self.dimension = 0 + for i in range(len(distributions)): + if isinstance(distributions[i], DistributionContinuous1D): + self.dimension += 1 + elif isinstance(distributions[i], JointIndependent): + self.dimension += len(distributions[i].marginals) + + # Initialize attributes + self.alpha: float = np.nan + """Normalized gradient vector in :math:`\\textbf{U}` space""" + self.alpha_record: list = [] + """Record of :math:`\\alpha=\\frac{\\nabla G(u)}{||\\nabla G(u)||}`""" + self.beta_record: list = [] + """Record of Hasofer-Lind reliability index that defines the feasibility criteria + :math:`||\\textbf{U}||=\\beta_{HL}`""" + self.design_point_u: list = [] + """Design point in the standard normal space :math:`\\textbf{U}`""" + self.design_point_x: list = [] + """Design point in the parameter space :math:`\\textbf{X}`""" + self.error_record: list = [] + """Record of the final error defined by + :math:`error_u = ||u_{new} - u||` and :math:`error_{\\nabla G(u)} = || \\nabla G(u_{new}) - \\nabla G(u)||`""" + self.iteration_record: list = [] + """Record of the number of iterations before algorithm termination""" + self.failure_probability_record: list = [] + """Record of the probability of failure defined by :math:`p_{fail} = \\Phi(-\\beta_{HL})`""" + self.state_function: list = [] + """State function :math:`G(u)` evaluated at each step in the optimization""" + + 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: + self.run(seed_x=self.seed_x) + + def run(self, seed_x: Union[list, np.ndarray] = None, seed_u: Union[list, np.ndarray] = None): + """Runs the inverse FORM algorithm. + + :param seed_x: Point in the parameter space :math:`\mathbf{X}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. + :param seed_u: Point in the uncorrelated standard normal space :math:`\mathbf{U}` to start from. + Only one of :code:`seed_x` or :code:`seed_u` may be provided. + If neither is provided, the zero vector in :math:`\mathbf{U}` space is the seed. + """ + self.logger.info('UQpy: Running InverseFORM...') + 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') + + # Allocate u and the gradient of G(u) as arrays + u = np.zeros([self.max_iterations + 1, self.dimension]) + state_function = np.full(self.max_iterations + 1, np.nan) + state_function_gradient = np.zeros([self.max_iterations + 1, self.dimension]) + + # Set initial seed. If both seed_x and seed_u are None, the initial seed is a vector of zeros in U space. + if seed_u is not None: + u[0, :] = seed_u + elif seed_x is not None: + self.nataf_object.run(samples_x=seed_x.reshape(1, -1), jacobian=False) + seed_z = self.nataf_object.samples_z + u[0, :] = Decorrelate(seed_z, self.nataf_object.corr_z).samples_u + + converged = False + iteration = 0 + while (not converged) and (iteration < self.max_iterations): + self.logger.info(f'Number of iteration: {iteration}') + if iteration == 0: + if seed_x is not None: + x = seed_x + else: + seed_z = Correlate(samples_u=u[0, :].reshape(1, -1), corr_z=self.nataf_object.corr_z).samples_z + self.nataf_object.run(samples_z=seed_z.reshape(1, -1), jacobian=True) + x = self.nataf_object.samples_x + else: + z = Correlate(u[iteration, :].reshape(1, -1), self.nataf_object.corr_z).samples_z + self.nataf_object.run(samples_z=z, jacobian=True) + x = self.nataf_object.samples_x + self.logger.info(f'Design Point U: {u[iteration, :]}\nDesign Point X: {x}\n') + state_function_gradient[iteration + 1, :], qoi, _ = self._derivatives(point_u=u[iteration, :], + point_x=x, + runmodel_object=self.runmodel_object, + nataf_object=self.nataf_object, + df_step=self.df_step, + order='first') + self.logger.info(f'State Function: {qoi}') + state_function[iteration + 1] = qoi + + alpha = state_function_gradient[iteration + 1] + alpha /= np.linalg.norm(state_function_gradient[iteration + 1]) + u[iteration + 1, :] = -alpha * self.beta + + error_u = np.linalg.norm(u[iteration + 1, :] - u[iteration, :]) + error_gradient = np.linalg.norm(state_function_gradient[iteration + 1, :] + - state_function_gradient[iteration, :]) + + converged_u = True if (self.tolerance_u is None) \ + else (error_u <= self.tolerance_u) + converged_gradient = True if (self.tolerance_gradient is None) \ + else (error_gradient <= self.tolerance_gradient) + converged = converged_u and converged_gradient + + if not converged: + iteration += 1 + + if iteration >= self.max_iterations: + self.logger.info(f'UQpy: Maximum number of iterations {self.max_iterations} reached before convergence') + self.alpha_record.append(alpha) + self.beta_record.append(self.beta) + self.design_point_u.append(u[iteration, :]) + self.design_point_x.append(x) + self.error_record.append((error_u, error_gradient)) + self.iteration_record.append(iteration) + self.failure_probability_record.append(self.p_fail) + self.state_function.append(state_function) diff --git a/src/UQpy/reliability/taylor_series/__init__.py b/src/UQpy/reliability/taylor_series/__init__.py index c5af5bab..e2797682 100644 --- a/src/UQpy/reliability/taylor_series/__init__.py +++ b/src/UQpy/reliability/taylor_series/__init__.py @@ -1,3 +1,4 @@ from UQpy.reliability.taylor_series.FORM import FORM from UQpy.reliability.taylor_series.SORM import SORM +from UQpy.reliability.taylor_series.InverseFORM import InverseFORM from UQpy.reliability.taylor_series.baseclass.TaylorSeries import TaylorSeries diff --git a/tests/unit_tests/reliability/example_7_2.py b/tests/unit_tests/reliability/example_7_2.py new file mode 100644 index 00000000..cbeb6f33 --- /dev/null +++ b/tests/unit_tests/reliability/example_7_2.py @@ -0,0 +1,18 @@ +import numpy as np + + +def performance_function(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 diff --git a/tests/unit_tests/reliability/test_inverse_form.py b/tests/unit_tests/reliability/test_inverse_form.py new file mode 100644 index 00000000..1c0e753e --- /dev/null +++ b/tests/unit_tests/reliability/test_inverse_form.py @@ -0,0 +1,105 @@ +import os +import pytest +import numpy as np +from scipy import stats +from UQpy.distributions import Normal +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.PythonModel import PythonModel +from UQpy.reliability.taylor_series import InverseFORM + + +@pytest.fixture +def inverse_form(): + """Example 7.2 from Chapter 7 of X. Du 2005 + + Distributions are :math:`P_X \\sim N(500, 100)`, :math:`P_Y \\sim N(1000, 100)` + Solution from the reference is :math:`u^*=(1.7367, 0.16376)`. + Tolerances of :math:`1e-5` are used to ensure convergence to level of precision given by Du. + """ + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + return InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=0.04054, + tolerance_u=1e-5, + tolerance_gradient=1e-5) + + +def test_no_seed(inverse_form): + inverse_form.run() + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_seed_x(inverse_form): + seed_x = np.array([625, 900]) + inverse_form.run(seed_x=seed_x) + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_seed_u(inverse_form): + seed_u = np.array([2.4, -1.0]) + inverse_form.run(seed_u=seed_u) + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-4) + + +def test_both_seeds(inverse_form): + """Expected behavior is to raise ValueError and inform user only one input may be provided""" + seed_x = np.array([1, 2]) + seed_u = np.array([3, 4]) + with pytest.raises(ValueError, match='UQpy: Only one input .* may be provided'): + inverse_form.run(seed_u=seed_u, seed_x=seed_x) + + +def test_neither_tolerance(): + """Expected behavior is to raise ValueError and inform user at least one tolerance must be provided""" + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + with pytest.raises(ValueError, match='UQpy: At least one tolerance .* must be provided'): + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=0.04054, + tolerance_u=None, + tolerance_gradient=None) + + +def test_beta(): + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=None, + beta=-stats.norm.ppf(0.04054)) + inverse_form.run() + assert np.allclose(inverse_form.design_point_u, np.array([1.7367, 0.16376]), atol=1e-3) + + +def test_no_beta_no_pfail(): + """Expected behavior is to raise ValueError and inform the user exactly one in put must be provided""" + path = os.path.abspath(os.path.dirname(__file__)) + os.chdir(path) + python_model = PythonModel(model_script='example_7_2.py', + model_object_name='performance_function', + delete_files=True) + runmodel_object = RunModel(model=python_model) + distributions = [Normal(loc=500, scale=100), Normal(loc=1_000, scale=100)] + with pytest.raises(ValueError, match='UQpy: Exactly one input .* must be provided'): + inverse_form = InverseFORM(distributions=distributions, + runmodel_object=runmodel_object, + p_fail=None, + beta=None) +