RELIABILITY: PARALLEL SYSTEM¶

In this example, the failure probability of a parallel system is computed using a plain Monte Carlo simulation (MCS) and the First-Order Reliability Method (FORM). The results are then compared.

Package imports¶

In [1]:
from uqpylab import sessions
import numpy as np
from scipy.stats import multivariate_normal as mvn

Start a remote UQCloud session¶

In [2]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (4655029c527342c98aecb7d14c894525) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [3]:
uq.rng(2,'twister');

COMPUTATIONAL MODEL¶

The parallel system consists of two two-dimensional Resistance-Stress (R-S) limit state functions that are defined as follows: $$ g_1(r, s) = 1.2 r - 0.9 s; \quad g_2(r, s) = r - s $$

Create two MODEL objects based on the limit state functions using string with vectorized operation:

In [4]:
Model1Opts = { 
    'Type': 'Model', 
    'mString': '1.2*X(:,1) - 0.9*X(:,2)',
    'isVectorized': 1
}
myLimitState1 = uq.createModel(Model1Opts)

Model2Opts = { 
    'Type': 'Model', 
    'mString': 'X(:,1) - X(:,2)',
    'isVectorized': 1
}
myLimitState2 = uq.createModel(Model2Opts)

The full parallel system limit state function can be calculated by just taking the maximum of the two limit state functions:

In [5]:
ModelFullOpts = { 
    'Type': 'Model', 
    'mString': 'max(X(:,1) - X(:,2), 1.2*X(:,1) - 0.9*X(:,2))',
    'isVectorized': 1
}
myLimitStateFull = uq.createModel(ModelFullOpts)

PROBABILISTIC INPUT MODEL¶

The probabilistic input model consists of two independent Gaussian random variables:

$$ R \sim \mathcal{N} (3, 0.3),\, S \sim \mathcal{N} (2, 0.4) $$

Specify the marginals of the two input random variables:

In [6]:
InputOpts = {
    "Marginals": [
        {"Name": "R",               # Resistance
         "Type": "Gaussian",
         "Moments": [3.0 , 0.3]
        },
        {"Name": "S",               # Stress
         "Type": "Gaussian",
         "Moments": [2.0 , 0.4]
        }
    ]
}

Create an INPUT object based on the specified marginals:

In [7]:
myInput = uq.createInput(InputOpts)

RELIABILITY ANALYSIS¶

Failure event is defined as $g(\mathbf{x}) \leq 0$. The failure probability is then defined as $P_f = P[g(\mathbf{x}) \leq 0]$.

Monte Carlo simulation (MCS)¶

Select the Reliability module and the Monte Carlo simulation (MCS) method:

In [8]:
MCSOpts = {
    "Type": "Reliability",
    "Method":"MCS"
}

MCS is performed on the full parallel system limit state function.

Select the function:

In [9]:
uq.selectModel(myLimitStateFull['Name'])

Specify the maximum sample size:

In [10]:
MCSOpts["Simulation"] = {
    "MaxSampleSize": 2e6
}

Run the Monte Carlo simulation:

In [11]:
myMCSAnalysis = uq.createAnalysis(MCSOpts)

Store the probality of failure in a variable:

In [12]:
Pf_MC = myMCSAnalysis['Results']['Pf']

First-order reliability method (FORM)¶

Select the Reliability module and the FORM method:

In [13]:
FORMOpts = {
    "Type": "Reliability",
    "Method":"FORM"
}

FORM is performed on each component of the system.

First, select the first component (i.e., the first limit state function):

In [14]:
uq.selectModel(myLimitState1['Name'])

Run the FORM analysis on the first limit state function:

In [15]:
myFORMAnalysis1 = uq.createAnalysis(FORMOpts)

Then, perform the FORM analysis on the second limit state function:

In [16]:
uq.selectModel(myLimitState2['Name'])
myFORMAnalysis2 = uq.createAnalysis(FORMOpts)

Retrieve the reliability index of each component:

In [17]:
betaHL1 = myFORMAnalysis1['Results']['BetaHL']
betaHL2 = myFORMAnalysis2['Results']['BetaHL']

Calculate the unit vectors ($\alpha$) in the direction of the design point for each system component:

In [18]:
alpha1 = np.array(myFORMAnalysis1['Results']['Ustar']) / betaHL1 
alpha2 = np.array(myFORMAnalysis2['Results']['Ustar']) / betaHL2

The failure probability of the parallel system is calculated by:

$P_f = \Phi_2(-\beta, 0, R),$

where $\Phi_2$ is a bivariate Gaussian distribution, with zero mean, unit variance, and correlation matrix $R$, evaluated at point $-\beta$.

The correlation matrix $R$ is defined as follows:

$R= \left[\begin{matrix} 1 & \alpha_1 \cdot \alpha_2 \\ \alpha_1 \cdot \alpha_2 & 1 \end{matrix} \right]$

Finally, the failure probability based on the FORM method is computed as follows:

In [19]:
B = np.array([betaHL1, betaHL2])
mu = np.zeros(2)
R = np.array([[1, alpha1@alpha2],
              [alpha1@alpha2, 1]])
dist = mvn(mean=mu, cov=R)
Pf_FORM = dist.cdf(-B)

For more details about the derivation of the formula see, for example, Chapter 15 of Engineering Design Reliability Handbook, CRC Press 2004.

RESULTS COMPARISON¶

Compare the failure probabilities obtained by the MCS and FORM:

In [20]:
print(f'Results:\n Pf[MC]   = {Pf_MC}\n Pf[FORM] = {Pf_FORM}')

Results:
 Pf[MC]   = 0.000203
 Pf[FORM] = 0.00020347600872246252

Terminate the remote UQCloud session:¶

In [21]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 4655029c527342c98aecb7d14c894525 terminated.
Out[21]:
True