Safety Factor Calibration#
Pystra includes the class Calibration which can be used to obtain calibrated load safety factors, \(\gamma\), resistance safety factor, \(\phi\), and load combination factors, \(\psi\). This example demonstrates running safety and combination factor calibration problems using Pystra.
Example#
This demonstration problem is adopted from Example 4 pp 190-1 from Sorensen, J.D. (2004), Notes in Structural Reliability Theory And Risk Analysis, with the following error-fixes:
The loads described as \(Q_1\) and \(Q_2\), are not the point-in-time loads, but the annual maxima distributions, from which the point-in-time loads are to be inferred.
There is an error in the textbook, and the wind load \(Q_2\) is not taken as occurring \(r_2 = 360\) times per year as stated, but instead \(r_2=2\) per year.
With these adaptations, the results of the calibration here match the results specified in Sorensen (2004).
Import Libraries#
[1]:
import pystra as ra
import numpy as np
Define the limit state function for calibration#
The LSF to be supplied in LoadCombination can be any valid Pystra LSF with one additional argument: a scalar design parameter for a random variable (generally, the resistance). Note that: 1. The design parameter must be specified as a Pystra constant with any arbitrary default value, such as 1.0. 2. Only a single design parameter for the resistance random variable can be calibrated at a time.
Keeping with the LoadCombination example, the design parameter for resistance, denoted \(z\), is added to the LSF:
[2]:
def lsf(z, R, G, Q1, Q2, cg):
return z*R - (cg*G + 0.6*Q1 + 0.3*Q2)
Define the Load and Resistance distributions#
Next follow the creation of load combination objects, as explained for the LoadCombination object.
Define the annual max distributions
[3]:
Q1max = ra.Gumbel("Q1", 1, 0.2) # Imposed Load
Q2max = ra.Gumbel("Q2", 1, 0.4) # Wind Load
Next specifcy the parameters of inferred point-in-time parent distributions
[4]:
Q1pit = ra.Gumbel("Q1", 0.89, 0.2) # Imposed Load
Q2pit = ra.Gumbel("Q2", 0.77, 0.4) # Wind Load
Q_dict = {'Q1': {'max': Q1max, 'pit': Q1pit},
'Q2': {'max': Q2max, 'pit': Q2pit}}
And define any constants to be passed through
[5]:
cg = ra.Constant("cg", 0.4)
z = ra.Constant("z", 1) # Design parameter for resistance with arbitrary default value
Finally, define any other random variables in the problem
[6]:
Rdist = ra.Lognormal("R", 1.0, 0.15) # Resistance
Gdist = ra.Normal("G", 1, 0.1) # Permanent Load (other load variable)
Instantiate LoadCombination object#
Class Calibration requires specification of a LoadCombination object to run reliability analyses for different load cases.
Setup load combination cases#
For this problem, we’re interested in investigating the reliabilities for two cases: 1. Q1_max: Q1 is maximum and Q2 is the point-in-time distribution 2. Q2_max: Q1 is the point-in-time distribution and Q2 is maximum
[7]:
loadcombinations = {'Q1_max':['Q1'], 'Q2_max':['Q2']}
[8]:
lc = ra.LoadCombination(lsf=lsf,
dict_dist_comb=Q_dict,
list_dist_resist=[Rdist],
list_dist_other=[Gdist],
list_const=[z, cg],
dict_comb_cases=loadcombinations)
Note that Pystra categorizes the variables for the calibration problem as follows: 1. Combination load variables: time variant load effects with a point-in-time distribution and maximum distribution 2. Resistance variables 3. Other load variables (optional): time invariant load effects 4. Constants: including the calibration parameter
Specify Nominal (or characteristic) values#
In calibration problems, the safety and combination factors are calibrated with respect to a given nominal value of load and resistance distributions. Generally, these nominal values correspond to the characteristic value of the load or resistance, e.g. 95-percentile or 5-percentile value.
Class Calibration requires specifying the nominal values as Python dictionary.
For this problem, the combination loads use a 98-percentile value for nominal, and a 5-percentile for resistance; the permanent load nominal value is considered at its mean. These are found using the ppf function in each Pystra distribution object:
[9]:
Qk = np.array([Q1max.ppf(0.98), Q2max.ppf(0.98)])
Gk = np.array([Gdist.ppf(0.5)])
Rk = np.array([Rdist.ppf(0.05)])
rvs_all = ['R', 'G', 'Q1', 'Q2', 'Q3']
dict_nom = dict(zip(rvs_all, np.concatenate([Rk, Gk, Qk])))
Specify target reliability index, \(\beta_T\)#
In calibration problems, the safety and combination factors are calibrated to obtain a lower bound safety indicated by a target reliablity index, \(\beta_T\).
Class Calibration requires specifying the target reliability index as a floating point integer.
[10]:
betaT = 4.3
This corresponds to a probability of failure of:
[11]:
ra.StdNormal.cdf(-betaT)
[11]:
8.539905471005582e-06
Calibrate safety and combination factors#
Class Calibration structures the calibration problem into three parts: 1. Calibration: this part involves projecting the design point of the LSF to correspond to the specified \(\beta_T\). 2. Estimation: this part involves estimating and calibrating the safety and combination factors, \(\phi\), \(\gamma\), and \(\psi\) using the projected design point. These factors correspond to \(\beta_T\) when utilized in structural design code equation,
\(z\phi R_n \geq \sum_i \psi_i \gamma_i S_{ni}\), where \(S_{ni}\) is the nominal value of the load effect \(S_i\). 3. Design check: this part involves utilizing the estimates of \(\phi\), \(\gamma\), and \(\psi\) for the given load effects to obtain the maximum required design resistance, i.e. \(z \cdot R_n\). Then the reliabilities (or safety) corresponding to the maximum required design resistance are checked for all load cases to ensure that they are greater
than \(\beta_T\). This design check is explained in more detail later in this tutorial.
Calibration algorithms#
Class Calibration implements two algorithms for calibration: optimize (default) and alpha.
Algorithm optimize (default)#
The optimize algorithm calibrates the resistance design parameter, \(z\), to obtain the design point on the LSF corresponding to \(\beta_T\) by utilizing SciPy’s inbuilt optimization algorithms.
Algorithm alpha#
The alpha algorithm calibrates the resistance design parameter, \(z\), by projecting the design point using the FORM \(\alpha\) estimates to obtain the design point on the LSF corresponding to \(\beta_T\).
Estimation algorithms#
Class Calibration implements two algorithms for estimation of safety and combination factors: coeff and matrix (default).
Algorithm matrix (default)#
The matrix algorithm utilizes the calibrated design points, \(X^*_{\beta_T}\), to estimate \(\phi\), \(\gamma\), and \(\psi\) by formulating a set of simultaneous equations using the structural design code equation, \(z\phi R_k \geq \sum_i \psi_i \gamma_i S_{ki}\), where \(S_{ki}\).
For more details, see Caprani, C. and Khan, M. S., Determination of load combination factors for the assessment of existing bridges, ICASP14 (July 2023).
The estimates of \(\phi\) and \(\gamma\) are obtained by dividing \(X^*_{\beta_T}\) with the corresponding nominal values as outlined earlier. To estimate \(\psi\), the structural design code equation is rearranged to obtain a system of simultaneous equation at the limiting condition,
Using numpy libraries to solve the system of simultaenous equations, estimates of \(\psi\) are obtained.
This matrix estimation algorithm provides unique values of \(\gamma\) and \(\psi\), while the \(\phi\) values can differ per load case. Calibration.run() takes an optional argument set_max (default value is False). If set_max=True, then the coeff algorithm sets the resulting dataframe of \(\phi\) to their maximum value per load case.
Algorithm coeff#
The coeff algorithm utilizes the calibrated design points, \(X^*_{\beta_T}\), to estimate \(\phi\), \(\gamma\), and \(\psi\) by comparing the coefficients. In this comparison, the design points values are mapped to the corresponding factors and nominal values, and the values of the corresponding factors are estimated, as demonstrated below. For more details refer to pp 136-42 of Sorensen, J.D. (2004), Notes in Structural Reliability Theory And Risk Analysis.
For the example problem, the design points recovered from the projection step for the two load cases are:
Dividing by the corresponding nominal values and considering only the random variables, we get,
By comparing the coefficients in the above matrix, estimates of \(\phi\), \(\gamma\), and \(\psi\) are obtained. This estimation algorithm provides unique values of \(\gamma\), while the \(\phi\) and \(\psi\) values can differ per load case. Calibration.run() takes an optional argument set_max (default value is False). If set_max=True, then the coeff algorithm sets the resulting dataframe of \(\psi\) or \(\phi\) to their maximum value per load
case.
Design Check#
Calibration implements methods to check whether the a design conducted on the basis of the estimated safety factors, actually achieves a minimum safety corresponding to \(\beta_T\) for all load cases. This check is implemented as follows. The calibrated estimates of \(\phi\), \(\gamma\), and \(\psi\) can be checked to ensure that a design based on them actually achieves a lower bound safety corresponding to \(\beta_T\) as follows: 1. Using \(\phi\), \(\gamma\),
\(\psi\), and the specified nominal values, estimate the design multiplier for resistance, \(z\) for each load case, \(j\) using the LSF. 2. Using the design multiplier, i.e. maximum \(z\), and corresponding load and resistance distributions, do forward reliability analysis to obtain the reliability index for each load case, \(\beta_j\). 3. Check \(\beta_j \geq \beta_T\).
Pystra implements two methods for the design check: calibration.get_design_param_factor() and calibration.calc_beta_design_param() for above step 1 and step 2, respectively.
Reliability Calibration (using coeff and optimize)#
Perform Calibration#
This configuration mirrors that of Sorensen’s (2004) solution noted above.
[12]:
calib = ra.Calibration(lc, target_beta=betaT, dict_nom_vals=dict_nom, calib_var='z',
est_method="coeff", calib_method="optimize",
print_output=False)
calib.run()
calib.print_detailed_output()
======================================================
X* =
R G Q1 Q2 z
Q1_max 0.66 1.04 1.62 2.02 3.04
Q2_max 0.66 1.04 1.51 2.25 3.05
phi =
R
Q1_max 0.85
Q2_max 0.85
gamma =
G Q1 Q2
Q1_max 1.04 1.07 1.1
Q2_max 1.04 1.07 1.1
psi =
G Q1 Q2
Q1_max 1.0 1.00 0.9
Q2_max 1.0 0.93 1.0
======================================================
The above results are in agreement with those presented in Example 4 in Sorensen, J.D. (2004), Notes in Structural Reliability Theory And Risk Analysis.
Design Check#
Here we confirm that the reliability target is met, as explained above.
[13]:
design_z = calib.get_design_param_factor()
design_beta = calib.calc_beta_design_param(np.max(design_z))
print(f"Design reliabilities = {design_beta.round(2)}")
print(f"Design Check = {design_beta.round(2)>=betaT}")
Design reliabilities = [4.31 4.3 ]
Design Check = [ True True]
Similarly, the calibration can be done using other algorithms
Reliability Calibration (using matrix and optimize)#
[14]:
calib = ra.Calibration(lc, target_beta=betaT, dict_nom_vals=dict_nom, calib_var='z',
est_method="matrix", calib_method="optimize")
calib.run()
calib.print_detailed_output()
======================================================
X* =
R G Q1 Q2 z
Q1_max 0.66 1.04 1.62 2.02 3.04
Q2_max 0.66 1.04 1.51 2.25 3.05
phi =
R
Q1_max 0.85
Q2_max 0.85
gamma =
G Q1 Q2
Q1_max 1.04 1.07 1.1
Q2_max 1.04 1.07 1.1
psi =
G Q1 Q2
Q1_max 1.0 1.00 0.9
Q2_max 1.0 0.93 1.0
======================================================
[15]:
design_z = calib.get_design_param_factor()
design_beta = calib.calc_beta_design_param(np.max(design_z))
print(f"Design reliabilities = {design_beta.round(2)}")
print(f"Design Check = {design_beta.round(2)>=betaT}")
Design reliabilities = [4.31 4.3 ]
Design Check = [ True True]
Reliability Calibration (using matrix and alpha)#
[16]:
calib = ra.Calibration(lc, target_beta=betaT, dict_nom_vals=dict_nom, calib_var='z',
est_method="matrix", calib_method="alpha")
calib.run()
calib.print_detailed_output()
======================================================
X* =
R G Q1 Q2 z
Q1_max 0.66 1.04 1.62 2.02 3.04
Q2_max 0.66 1.04 1.51 2.25 3.05
phi =
R
Q1_max 0.85
Q2_max 0.85
gamma =
G Q1 Q2
Q1_max 1.04 1.07 1.1
Q2_max 1.04 1.07 1.1
psi =
G Q1 Q2
Q1_max 1.0 1.00 0.9
Q2_max 1.0 0.93 1.0
======================================================
[17]:
design_z = calib.get_design_param_factor()
design_beta = calib.calc_beta_design_param(np.max(design_z))
print(f"Design reliabilities = {design_beta.round(2)}")
print(f"Design Check = {design_beta.round(2)>=betaT}")
Design reliabilities = [4.31 4.3 ]
Design Check = [ True True]
Safety Factor Calibration for Non Linear Limit State Functions (LSFs)#
This tutorial demonstrates how Pystra can calibrate safety and combination factors for Non-Linear limit state functions.
Import Libraries#
[18]:
import pystra as ra
import numpy as np
Define the limit state function for calibration#
The LSF to be supplied in LoadCombination can be any valid Pystra LSF with one additional argument: a scalar design parameter for a random variable (generally, the resistance). Note that: 1. The design parameter must be specified as a Pystra constant with any arbitrary default value, such as 1.0. 2. Only a single design parameter for the resistance random variable can be calibrated at a time.
Keeping with the LoadCombination example, the design parameter for resistance, denoted \(z\), is added to the LSF. The LSF contains additional random variables wR and wS to account for the uncertainties within the resistance and the load effects, respectively. The non-linearity of the LSF is evident as the degree of random variables in each term of the LSF is two.
[19]:
def lsf(z, wR, wS, R, Q1, Q2):
gX = z * wR * R - wS * (Q1 + Q2)
return gX
Define the Load and Resistance distributions#
Next follow the creation of load combination object, as explained for the LoadCombination class.
[20]:
wR = ra.Lognormal("wR", 1.0, 0.05)
wS = ra.Lognormal("wS", 1.0, 0.10)
R = ra.Normal("R", 60, 6) # [units]
Q1_max = ra.Normal("Q1", 30, 3) # [units]
Q2_max = ra.Normal("Q2", 20, 2) # [units]
Q1_pit = ra.Normal("Q1", 15, 3) # [units]
Q2_pit = ra.Normal("Q2", 10, 2) # [units]
z = ra.Constant("z", 1)
Set up Load Combinations#
Specify nominal values & combinations#
[21]:
rvs_all = ['wR', 'wS', 'R', 'Q1', 'Q2']
dict_nom = dict(zip(rvs_all, np.array([1.0, 1.0, R.ppf(0.05),
Q1_max.ppf(0.95),
Q2_max.ppf(0.95)])))
Q_dict = {'Q1': {'max': Q1_max, 'pit': Q1_pit},
'Q2': {'max': Q2_max, 'pit': Q2_pit}}
[22]:
loadcombinations = {'Q1_max':['Q1'], 'Q2_max':['Q2']}
Instantiate LoadCombination object#
[23]:
lc = ra.LoadCombination(lsf,
dict_dist_comb=Q_dict,
list_dist_other=[wS],
list_dist_resist=[R, wR],
list_const = [z],
dict_comb_cases=loadcombinations)
Note that since wR is associated with the resistance, it is specified in list_dist_resist; while wS is specified in list_dist_other as it is associated with the load effects but it’s not a load combination variable.
Safety Factor calibration#
[24]:
betaT = 3.7
calib1 = ra.Calibration(lc, target_beta=betaT, dict_nom_vals=
dict_nom, calib_var='z',
est_method="matrix", calib_method="optimize")
calib1.run()
calib1.print_detailed_output()
======================================================
X* =
R wR wS Q1 Q2 z
Q1_max 44.40 0.95 1.2 33.81 11.69 1.30
Q2_max 44.76 0.95 1.2 19.16 21.85 1.16
phi =
R wR
Q1_max 0.89 0.95
Q2_max 0.89 0.95
gamma =
wS Q1 Q2
Q1_max 1.2 0.97 0.94
Q2_max 1.2 0.97 0.94
psi =
wS Q1 Q2
Q1_max 1.0 1.00 0.54
Q2_max 1.0 0.57 1.00
======================================================
Design Check#
[25]:
design_z1 = calib1.get_design_param_factor()
design_beta1 = calib1.calc_beta_design_param(np.max(design_z1))
print(f"Design reliabilities = {design_beta1.round(2)}")
print(f"Design Check = {design_beta1.round(2)>=betaT}")
Design reliabilities = [3.7 4.29]
Design Check = [ True True]