Introductory Tutorial#
This is a short introduction on how to use pystra. The tutorial above is also available on GitHub under example.py.
An example reliability model#
Consider the following random variables:
Additionally those variables are related to each other. Therefore the correlation matrix \({\bf C}\) is given:
Now, we like to compute the reliability index \(\beta\) and the failure probability \(P_f\), by given limit state function \(g(\gamma, X_1,X_2,X_3)\):
where \(\gamma\) is a real constant. For this example, let \(\gamma = 1\)
Establish the model#
Before we start with the modeling, we have to import the pystra package and other relevant packages.
[1]:
import pystra as ra
import numpy as np
Two ways to define the limit state function are available: + Direct in the main code, + as a separate function.
In the first case the input will look like:
[2]:
# Define limit state function
# - case 1: define directly
limit_state = ra.LimitState(lambda g,X1,X2,X3: g - X2*(1000*X3)**(-1) - (X1*(200*X3)**(-1))**2)
and in the second case like this:
[3]:
# Define limit state function
# - case 2: use predefined function
def example_limitstatefunction(g, X1, X2, X3):
"""
example limit state function
"""
return g - X2 * (1000 * X3) ** (-1) - (X1 * (200 * X3) ** (-1)) ** 2
limit_state = ra.LimitState(example_limitstatefunction)
Notice, here the function example_limitstatefunction has be defined in advance as a separate function. This case can be useful if the limit state function is quiet complex or need more then one line to define it.
In the next step the stochastic model has to be initialized
[4]:
stochastic_model = ra.StochasticModel()
and the random variables have to be assigned. To define the random variables from (1) we can use following syntax:
[5]:
# Define random variables
stochastic_model.addVariable(ra.Lognormal("X1", 500, 100))
stochastic_model.addVariable(ra.Normal("X2", 2000, 400))
stochastic_model.addVariable(ra.Uniform("X3", 5, 0.5))
The first parameter is the name of the random variable. The name has to be a string and match the arguments in the limit state function, so the input looks like "X3".
By default, the next to values are the first and second moment of the distribution, here mean and standard deviation. If mean and standard deviation unknown but the distribution parameter known, then the input_type has to be changed.
For example random variable \(X_3\) is uniform distributed. Above we assume that \(X_3\) is defined by mean and standard deviation. But we can describe the distribution with the parameter \(a\) and \(b\). In this case the code will look like:
[6]:
X3 = ra.Uniform('X3',4.133974596215562, 5.866025403784438, 1)
to get the same results as before. To see which parameters are needed and in which order the must insert, refer to the Distributions API.
If the nominal value, bias, and coefficient of variation are instead known, then the random variable can be instantiated following this example:
[7]:
X2 = ra.Normal('X2',*500*1.00*np.array([1, 0.2]))
where nominal value is 500, bias is 1.00, and coefficient of variation is 0.2. Notice the initial * character is used to dereference the output array.
We will also define our constant using Constant:
[8]:
# Define constants
stochastic_model.addVariable( ra.Constant('g',1) )
To add the correlation matrix to our model:
[9]:
# Define Correlation Matrix
stochastic_model.setCorrelation( ra.CorrelationMatrix([[1.0, 0.3, 0.2],
[0.3, 1.0, 0.2],
[0.2, 0.2, 1.0]]) )
If the variables uncorrelated, you don’t have to add a correlation matrix to the model.
At this stage our model is complete defined and we can start the analysis.
Reliability Analysis#
To change some options, a object must be initialized which stores the customized options.
[10]:
options = ra.AnalysisOptions()
options.setPrintOutput(False)
To store the results from the analysis an object must be initialized ### FORM Analysis Now the code can be compiled and the FORM analysis will be preformed. In this example we will get following results:
[11]:
# initialize analysis obejct
Analysis = ra.Form(
analysis_options=options,
stochastic_model=stochastic_model,
limit_state=limit_state,
)
Analysis.run() # run analysis
If we don’t like to see the results in the terminal the option setPrintOutput(False) has set to be False. There are also some other options which can be modified.
To use the results for further calculations, plots etc. the results can get by some getter methods
[12]:
# Some single results:
beta = Analysis.getBeta()
failure = Analysis.getFailure()
There is also the possibility to output more detailed results using showDetailedOutput():
[13]:
Analysis.showDetailedOutput()
======================================================
FORM
======================================================
Pf 3.9717297753e-02
BetaHL 1.7539761407
Model Evaluations 164
------------------------------------------------------
Variable U_star X_star alpha
X1 1.278045 631.504135 +0.728414
X2 0.407819 2310.352495 +0.232354
X3 -1.129920 4.517374 -0.644534
g --- 1.000000 ---
======================================================
SORM Analysis#
SORM Analysis#
A Second-Order Reliability Method (SORM) can also be performed, passing in the results of a FORM analysis object if it exists. For efficiency, we can pass the FORM results object if it exists, otherwise it will be called automatically.
[14]:
sorm = ra.Sorm(
analysis_options=options,
stochastic_model=stochastic_model,
limit_state=limit_state,
form=Analysis,
)
sorm.run()
Similar to FORM, we can also get more detailed output for diagnostics:
[15]:
sorm.showDetailedOutput()
======================================================
FORM/SORM
======================================================
Pf FORM 3.9717297753e-02
Pf SORM Breitung 3.2229053013e-02
Pf SORM Breitung HR 3.1158626135e-02
Beta_HL 1.7539761407
Beta_G Breitung 1.8489979687
Beta_G Breitung HR 1.8640317038
Model Evaluations 180
------------------------------------------------------
Curvature 1: -0.04143130874014485
Curvature 2: 0.36356407428350895
------------------------------------------------------
Variable U_star X_star alpha
X1 1.278045 631.504135 +0.728414
X2 0.407819 2310.352495 +0.232354
X3 -1.129920 4.517374 -0.644534
g --- 1.000000 ---
======================================================
in which HR refers to the Hohenbichler-Rackwitz modification to Breitung’s formula.
Distribution Analysis#
[16]:
da = ra.DistributionAnalysis(
analysis_options=options,
stochastic_model=stochastic_model,
limit_state=limit_state,
)
da.run()
Crude Monte Carlo Simulation#
[17]:
cmc = ra.CrudeMonteCarlo(
analysis_options=options,
stochastic_model=stochastic_model,
limit_state=limit_state,
)
cmc.run()
Importance Sampling#
[18]:
ismc = ra.ImportanceSampling(
analysis_options=options,
stochastic_model=stochastic_model,
limit_state=limit_state,
)
ismc.run()
Results#
[19]:
beta = Analysis.getBeta()
failure = Analysis.getFailure()
print(f"Beta is {beta}, corresponding to a failure probability of {failure}")
Beta is 1.7539761407409655, corresponding to a failure probability of [0.0397173]