OpenSees Cloud

OpenSees AMI

The Toy Reliability Problem

Original Post - 23 Jun 2024 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.


In Fundamentals of Structural Mechanics, Hjelmstad uses the “little BVP” to introduce boundary value problems. Chopra uses a two-story rigid shear frame throughout Dynamics of Structures. And in Structural and System Reliability, Der Kiureghian uses a “toy problem” for various types of reliability analyses.

The toy reliability problem is defined by two random variables with the following limit state function

\[g = X_1^2 - 2X_2\]

where the mean and standard deviation of X1 are 10 and 2, respectively, and 20 and 5 for X2. The random variables have correlation coefficient 0.5.

This post will focus on using OpenSeesPy to perform first order reliability method (FORM) analysis of the toy limit state function where the random variables are Nataf-distributed with X1 lognormal and X2 Type I largest value (Gumbel). This analysis is Example 6.2 from ADK’s textbook.

On a side note, some of the textbook’s examples and content were developed around the time Terje Haukaas implemented reliability modules for OpenSees. Reading the book brought back some memories! Full details of the OpenSees reliability implementation, based on Tcl, can be found in this PEER report. But now, thanks to Minjie, the reliability commands are available in OpenSeesPy.

Defining the random variables for the toy problem is straightforward using the specified means, standard deviations, and correlation.

ops.wipe()
ops.wipeReliability()
 
ops.randomVariable(1,'lognormal','-mean',10.0,'-stdv',2.0)
ops.randomVariable(2,'gumbel','-mean',20.0,'-stdv',5.0)
 
ops.correlate(1,2,0.5)

Although tags of 1 and 2 are used here, the random variable tags do not have to be sequential.

Next, we define the transformation of random variable realizations between the standard and normal spaces. Whether the random variables are jointly normal or Nataf-distributed, the Nataf probability transformation will work.

ops.probabilityTransformation('Nataf')

At the moment, all random variables must be mapped to parameter objects in OpenSeesPy, which is a little different from reliability analysis in Tcl-based OpenSees. For this toy problem, we can define two individual parameters.

#         paramTag    type     rvTag
ops.parameter(1,'randomVariable',1)
ops.parameter(2,'randomVariable',2)

The parameter tags do not have to be the same as the random variable tags.

When there are many random variables, defining parameters one-by-one can become tedious. Instead, we can loop over the random variables and define parameters automatically.

for rv in ops.getRVTags():
    ops.parameter(rv,'randomVariable',rv)

Next, we can define the limit state function, also known as a performance function. The par variable is a special dictionary defined internally within the OpenSeesPy interpreter. The keys to par are the random variable tags and the values are the current realizations of the random variables. We define the limit state function \(g = X_1^2 - 2X_2\) in a string using par, then indicate we will use Python to evaluate the string.

ops.performanceFunction(1, "par[1]**2-2*par[2]")
 
ops.functionEvaluator('Python')

Every time the limit state function is evaluated internally during an OpenSees reliability analysis, the par dictionary will contain the current values of the random variables.

Also note that the limit state function tags do not have to be sequential.

Next, we define how the gradient of the limit state function will be evaluated during the search for the design point. The first option is FiniteDifference, which uses parameter perturbations of 1/1000 by default.

ops.gradientEvaluator('FiniteDifference')

The second option is the Implicit gradient evaluator, which allows the user to define “exact” expressions for the derivative of the limit state function with respect to the random variables. For this toy problem, the gradients are \(\partial g/\partial X_1 = 2X_1\) and \(\partial g/\partial X_2 = -2\).

ops.gradientEvaluator('Implicit')
 
#                         lsfTag rvTag expression
ops.gradPerformanceFunction(1,1,"2*par[1]")
ops.gradPerformanceFunction(1,2,"-2")

While the difference between the FiniteDifference and Implicit gradient evaluators is inconsequential for this toy problem, the choice can be quite consequential for finite element reliability analysis because the Implicit gradient evaluator can utilize the direct differentiation method (DDM) in the gradient calculations.

With the random variables and the limit state function and its gradients defined, we need to specify how the design point will be found. The reliability analysis options are much like the regular finite element analysis options in OpenSees–you get to pick and choose.

First is the starting point, the step size, and the search algorithm. The mean values of the random variables is usually the best starting point. You can also start from all random variables equal zero, but that’s usually pretty far from the design point. For this toy problem, fixed step size is OK, but for other problems, you may want to use the Armijo rule.

ops.startPoint('Mean')
ops.stepSizeRule('Fixed')

Next is the search algorithm and the convergence check. The improved Hasofer-Lind and Rackwitz-Fiessler (iHLRF) algorithm is the most common iterative search algorithm, but there are others implemented in OpenSees, such as Polak-He, gradient projection, and sequential quadratic programming (SQP). The standard reliability convergence check tests the value of the limit state function relative to its initial value (-e1) and the collinearity of the gradient vector and the vector of standard normal random variables (-e2). Both convergence criteria must be satisfied.

ops.searchDirection('iHLRF')
ops.reliabilityConvergenceCheck('Standard','-e1',1e-3,'-e2',1e-3,'-print',1)

With a search algorithm and convergence check, we can tell how to find the design point. The search with step size and direction (StepSearch) is the only option. Here we specify the maximum number of iterations as well.

ops.findDesignPoint('StepSearch','-maxNumIter',10)

Finally, we perform the FORM analysis and give an output file name.

ops.runFORMAnalysis('toy-reliability.out')

Running the analysis gives the output shown below. The analysis takes six steps to converge to the design point where the reliability index is 2.69 with 0.362% as the corresponding probability of failure.

Reliability analysis output

In addition to screen output, the output file generated by the FORM analysis contains the number of steps required to find the design point as well as the reliability index, probability of failure, and the value of the random variables at the design point.

Reliability analysis output file

Extracting information from the CalREL-esque formatted output file can be cumbersome. Instead, you can use special dictionaries to work with the results directly in your Python scripts.

The code below reports all the essential information for the toy reliability problem.

for lsf in ops.getLSFTags():
    print(f'Limit State Function: {lsf}')
    print(f'beta = {ops.betaFORM[lsf]}')
    print(f'pf = {ops.pfFORM[lsf]}')
    for rv in ops.getRVTags():
         print(f'RV tag = {rv}')
         print(f'  x* = {ops.designPointXFORM[lsf,rv]}')
         print(f'  u* = {ops.designPointUFORM[lsf,rv]}')
         print(f'  alpha = {ops.alphaFORM[lsf,rv]}')
         print(f'  gamma = {ops.gammaFORM[lsf,rv]}')

Custom output for reliability analysis

Of course, with Python, you’ll want to use numpy, pandas, and other libraries to do much more fancy things than print the reliability analysis results to the screen!


© 2019-2025 by Michael H. Scott. All rights reserved.