---
title: "TesiproV – Reliability Analysis Methods for Structural Engineering"
author: "Til Lux, Konstantin Nille-Hauf, Tania Feiri, Marcus Ricker"
date: "`r format(Sys.Date(), '%Y-%m-%d')`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{TesiproV – Reliability Analysis Methods for Structural Engineering}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This document describes how the TesiproV package can be used. The package was developed as part of the "TesiproV" research project (sponsored by BMWi, the German Federal Ministry for Economic Affairs and Energy) and is used for the probabilistic safety assessment of structural components and systems. The source code of the package is freely available as open source.


# Introduction

The **TesiproV** package provides a flexible and extensible framework for *probabilistic reliability analysis* of structural components and systems.  
It enables engineers and researchers to quantify the safety level of structures by evaluating the probability of failure under uncertain loads, material properties and geometric parameters.

In classical deterministic design procedures, reliability is ensured through partial safety factors that implicitly account for uncertainties.  
However, such approaches often lack transparency regarding the actual reliability achieved.  
Probabilistic methods overcome this limitation by explicitly modelling all relevant sources of uncertainty — e.g., variability in material strength, load effects or model uncertainties — as random variables with defined statistical distributions.

The fundamental concept behind reliability analysis is the *limit‑state function* $g(X) = 0$, which separates safe from failed states:

$$
g(X) = R(X) - E(X) = 0
$$

where $R(X)$ represents the resistance and $E(X)$ denotes the load effect.  
Failure occurs when $g(X) < 0$. The random vector $X = [X_1, X_2, …]$ contains all basic variables describing these uncertainties.

The probability of failure can then be expressed as:

$$
P_f = P[g(X) < 0]
$$

and converted into a *reliability index* β via

$$
β = -Φ^{-1}(P_f)
$$

where Φ denotes the cumulative standard normal distribution.  
Both the reliability index and the probability of failure are measures used to assess and compare reliability.
The calculated probability of failure differs from the actual probability of failure, as the governing factor for the failure of a component or structure — human error — is not captured in the analysis.

TesiproV implements several algorithms to estimate these quantities efficiently:

- **MVFOSM** — Mean Value First Order Second Moment method; fast linear approximation based on moment statistics.
- **FORM** — First Order Reliability Method; iterative search for design point on failure surface.
- **SORM** — Second Order Reliability Method; curvature correction around design point.
- **Monte Carlo simulations** — direct numerical evaluation using random sampling or importance/subset sampling schemes.

These methods allow both single limit‑state analyses and complex system evaluations involving multiple interacting components.  
In addition, TesiproV supports parametric studies where selected input parameters are systematically varied to investigate sensitivity or perform probabilistic optimisation.

Independent of the solution algorithms, the effect of quality control on the distribution of a basic variable can be taken into account. Sampling plans with their operating characteristic (OC) curves are defined by `PROB_QUALCON`, and the posterior distribution of the controlled variable is derived by Bayesian updating. The updated variable is then used like any other basic variable in the reliability analysis.

Since version 0.9.5 the entire code base has been refactored to improve consistency between probabilistic objects (`PROB_BASEVAR`, `SYS_LSF`, `SYS_PROB`, …), error handling in limit‑state functions (`$check()`), transformations between mean/standard deviation and native distribution parameters (`$prepare()`), and parallel execution via the **future** package was added.

Version 1.0.0 is a consolidated release of the package, based on more than two years of use. It solves single limit-state functions as well as systems of several limit-state functions. Since this version, the influence of quality controls on the distribution of basic variables can also be taken into account (see the section on `PROB_QUALCON`).

Beyond its computational capabilities, TesiproV aims to provide a transparent educational tool for understanding reliability theory in structural engineering practice.  
Its object‑oriented design makes it straightforward to define new limit‑state functions, integrate additional distributions or link external solvers, e.g., for advanced Monte Carlo techniques.

The project was originally initiated at Hochschule Biberach in cooperation with the Institute of Structural Concrete at RWTH Aachen University within a BMWi research program (*Allgemeingültiges Verfahren zur Herleitung von Teilsicherheitsbeiwerten im Massivbau auf Basis probabilistischer Verfahren anhand ausgewählter Versagensarten (TesiproV) – Erstellung eines Richtlinienentwurfs*, FKZ 03TNK003E) and is now maintained by the Chair of Structural Concrete at TU Dortmund University.

Source code repository: <https://gitlab.tu-dortmund.de/ls-massivbau/tesiprov>  
License: [MIT](https://opensource.org/licenses/MIT)


# Installation

Official releases of the R package TesiproV are available on CRAN (<https://CRAN.R-project.org/package=TesiproV>) and can be installed in R as: 

```{r, eval = FALSE}
install.packages("TesiproV")
```

You can install a pre-release version of TesiproV directly from the GitLab server at TU Dortmund University:

```{r, eval = FALSE}
remotes::install_gitlab("ls-massivbau/tesiprov", host = "gitlab.tu-dortmund.de")
```

If you need any permissions, please contact <tesiprov.lmb.ab@tu-dortmund.de>.

Make sure to initialize the TesiproV package before using it:

```{r}
library("TesiproV")
```

## Conceptual background

Reliability analysis quantifies how likely a structure will perform safely throughout its lifetime despite inherent uncertainties.
Typical sources of uncertainty include variations in geometry (e.g., effective depth), material properties (e.g., concrete strength),
loading conditions (permanent vs variable actions), environmental influences and modelling simplifications.
Within TesiproV every uncertain quantity is represented by an object-oriented random variable definition.
Transformation between mean value / standard deviation and native distribution parameters is handled automatically,
ensuring consistent probabilistic input even when empirical data are incomplete.
System reliability can be evaluated either analytically using bounds (“SimpleBounds”) or numerically via simulation-based approaches.
Parallelisation through the parallel (base R) and future backend allows efficient computation even for large-scale Monte-Carlo analyses on multi-core machines.


# Basic building blocks

## Random variable objects

A basic random variable describes an uncertain quantity used in a limit-state function.
It is defined by its statistical distribution type (including parameters) or its statistical moments (mean value (*Mean*) and standard deviation (*Sd*)).
Supported distributions include:

| Distribution | Identifier | Package | Description |
|---------------|-------------|----------|--------------|
| Normal | `"norm"` | stats | Gaussian distribution defined by mean and standard deviation |
| Lognormal | `"lnorm"` | stats | Distribution of a variable whose logarithm is normally distributed |
| Shifted lognormal | `"slnorm"` | TesiproV | Lognormal distribution shifted by the lower bound x0 |
| Gumbel (Type I extreme value) | `"gumbel"` | evd | Used for modelling maxima or minima (extreme values) |
| Gamma | `"gamma"` | stats | Positive-valued model defined by shape and scale parameters |
| Exponential | `"exp"` | stats | Special case of gamma with shape = 1 |
| Beta | `"beta"` | stats | Bounded on [0,1], often used for proportions or probabilities |
| Weibull | `"weibull"` | stats | Flexible positive-valued model for material strength or lifetime data |
| Student-t | `"st"` | TesiproV | Heavy-tailed, marginal distribution of a normal variable with uncertain mean and standard deviation |
| Log-Student-t | `"lt"` | TesiproV | Heavy-tailed, marginal distribution of a lognormal variable with uncertain parameters |
| Empirical | `"emp"` | EnvStats | Empirical distribution based on given data |
| Binomial | `"binom"` | stats | Discrete count data, placeholder type |
| Uniform | `"unif"` | stats | Constant probability across an interval, placeholder type |


Additional placeholders exist for future extensions.

### Example – PROB_BASEVAR

```{r}
var_fck <- PROB_BASEVAR(
  Id = 1,
  Name = "f_ck",
  Description = "Concrete compressive strength",
  DistributionType = "lnorm",
  Mean = 47.35,
  Sd = 5.86
)
```


* Id - numeric, for identification purpose 
* Name - Name of the var, must be the same as in the function input
* Description - Description for the output protocol
* Package - Packagename, which provides the distribution function (dname,pname,qname,rname; name=DistributionType) 
* DistributionType - Name of the distributionfunction (without the d, p, q or r)
* Mean - Expected value (first statistical moment)
* Sd - standard deviation (square root of second statistical moment)
* Cov - coefficient of variation: Cov = Sd/Mean
* x0 - location parameter, used for shifted log-normal distribution "slnorm"

Internally TesiproV automatically derives all missing distribution parameters such as shape or scale parameters from the given expected value (mean) and standard deviation (square root of variance) if possible.

Basic variable using the example of concrete compressive strength.
As distribution function a logStudentT-distribution adjusted within the package is used. For this distribution a d-, p-, q- and rlt-function is available, therefore Package = "TesiproV".

```{r}
var.f_ck <- PROB_BASEVAR(Id = 1, Name = "f_ck", Description = "Concrete compressive strength", Package = "TesiproV", DistributionType = "lt", DistributionParameters = c(3.85, 0.09, 3, 10))
```

## Parametric random variables

To perform parameter studies over varying input values (e.g., different mean strengths or load levels), use `PARAM_BASEVAR`. It extends PROB_BASEVAR with fields defining which property should vary (`ParamType`) and which values are tested (`ParamValues`).

### Example – PARAM_BASEVAR

```{r}
pvar_d <- PARAM_BASEVAR(
  Id = 1,
  Name = "d",
  Description = "Effective depth",
  DistributionType = "norm",
  Sd = 10,
  ParamType = "Mean",
  ParamValues = c(seq(200, 800, 50), seq(1000, 3000, 100)) + 10
)
```

## Deterministic variables

If certain quantities are treated as deterministic constants rather than random variables, they can be modelled using `PROB_DETVAR`. Internally these are represented as normal distributions with infinitesimal standard deviation so that all algorithms remain compatible.

### Example PROB_DETVAR

```{r}
var.gamma_c <- PROB_DETVAR(Id = 1, Name = "gamma_c", Description = "PSF Concrete", Value = 1.5)
```

During each iteration of a parametric study the corresponding value from `ParamValues` is assigned automatically.


## Parametric deterministic variables

Frequently, deterministic variables in a parameter calculation also depend on the varying parameter size. In these cases, a parametric deterministic variable can be created. With each iteration of the parametric base variable, the next step of the parametric deterministic variable is also used.

In this case, a normal distribution with infinitesimal standard deviation is also used internally. 

Note: Parametric variables are only considered for the parametric system object (SYS_PARAM) while the normal system object (SYS_PROB) only uses the default value.

### Example PARAM_DETVAR

```{r}
var.V_Ed <- PARAM_DETVAR(Id = 1, Name = "V_Ed", Description = "applied shear force", ParamValues = c(1, 2, 3, 4, seq(5, 10, 0.5)))
```


## Quality control and Bayesian updating (PROB_QUALCON)

Quality controls (e.g. the conformity control of concrete or the inspection of dimensions on site) accept or reject a lot. An accepted lot is on average better than the prior assumption. TesiproV quantifies this effect: the object `PROB_QUALCON` describes the sampling plan and computes its OC curve, and `PROB_BASEVAR$updateBasevar()` derives the posterior distribution of the controlled basic variable.

### Principle

A lot is characterised by its defect fraction $\theta$, i.e. the share of units that violate a limit value. The OC curve gives the probability that the lot is accepted by the sampling plan:

$$
P_a(\theta)
$$

The distribution of the controlled variable is described by a Student-t (`"st"`) or log-Student-t (`"lt"`) distribution with the hyperparameters $c(m', s', n', v')$, which also express the uncertainty of the distribution parameters. Every parameter combination belongs to one defect fraction $\theta$ and is therefore weighted with its acceptance probability $P_a(\theta)$. The result is the outgoing distribution of the accepted lots. It is fitted by a new `"st"` or `"lt"` distribution. If the quality of the fit is not sufficient (`min_r2` in `updateBasevar()`), an empirical distribution is returned instead.

The prior variable must therefore be defined with `DistributionType = "st"` or `"lt"` and the four hyperparameters in `DistributionParameters`.

### Sampling plans

The plan is chosen with the field `ConformityType`:

* `"fck_EN-206"` with `ConformityParam = c(n)`: variable plan, conformity criterion of EN 206 for the compressive strength. A lot is accepted if the mean of `n` samples is at least $f_{ck} + 1.48 s$ and the smallest value at least $f_{ck} - 4$. The OC curve of this plan is simulated. The Monte Carlo sample size per defect fraction is set with the field `NSim` (default 100000).
* `"dim_k||n"` with `ConformityParam = c(n, c_accept)`: attribute plan with a binomial model. A lot is accepted if at most `c_accept` of `n` units are defective. It can be used for any characteristic that is checked against tolerance values, e.g. concrete cover or dimensions.

The fields `DefectLimitLower` and `DefectLimitUpper` only define when a single unit counts as defective. They are not the acceptance criterion, which is determined by the sampling plan alone.

* For `"fck_EN-206"` the characteristic strength $f_{ck}$ is given as `DefectLimitLower` (in the unit of the variable). `DefectLimitUpper` must stay `Inf`.
* For `"dim_k||n"` at least one of the two limits must be finite. Both limits are freely selectable and do not influence the OC curve, only the update.

Fields and methods:

* `var` - the `PROB_BASEVAR` the plan applies to.
* `NSim` - Monte Carlo sample size per defect fraction for the simulation of the OC curve (only for `"fck_EN-206"`, default 100000).
* `theta_grid` - defect fractions at which the OC curve is evaluated (default grid from 0.001 to 0.9).
* `Pa_grid` - acceptance probabilities for `theta_grid`, filled by `$getOC()`.
* `$getOC(force = FALSE)` - computes the OC curve. An existing curve is only recomputed with `force = TRUE`.
* `$plot()` - plots the OC curve.
* `$isStale()` - `TRUE` if the distribution of `var` has changed since `$getOC()` was run (the curve is then outdated), `NA` if no snapshot exists.
* `$summary()` - short printed overview.

If the OC curve of a plan is known from another source, it can be specified manually by passing `theta_grid` and `Pa_grid` of the same length. `ConformityParam` is then not needed. For `"fck_EN-206"` `DefectLimitLower` must still be set.

### Example - variable plan for the concrete strength

The compressive strength is modelled by a log-Student-t distribution with the hyperparameters $c(3.75, 0.105, 3, 10)$. The control follows EN 206 with 15 samples for $f_{ck} = 30$ N/mm².

```{r, message=FALSE}
set.seed(1)
fc <- PROB_BASEVAR(
  Id = 1, Name = "fc", Description = "Concrete compressive strength",
  DistributionType = "lt", DistributionParameters = c(3.75, 0.105, 3, 10)
)

plan <- PROB_QUALCON(
  Name = "EN206_n15", ConformityType = "fck_EN-206",
  ConformityParam = 15, NSim = 20000, DefectLimitLower = 30, var = fc
)
plan$getOC()
```

```{r, fig.width=6, fig.height=4}
plan$plot()
```

`updateBasevar()` returns a new `PROB_BASEVAR` object with the posterior distribution. The prior object stays unchanged. The message printed by `updateBasevar()` compares the mean value and the standard deviation of prior and posterior. For hyperparameter-defined variables the prior moments are exact for `"st"` with $v' > 2$. For `"lt"` they are not defined analytically and are therefore approximated numerically, which is marked in the message.

```{r}
fc_post <- fc$updateBasevar(plan)
```

The moments of the posterior are also available as fields of the new object:

```{r}
c(Mean = fc_post$Mean, Sd = fc_post$Sd)
```

Prior and posterior density can be compared with the argument `prior`:

```{r, fig.width=6, fig.height=4}
fc_post$plot(prior = TRUE)
```

### Example - attribute plan for the concrete cover

At most one of ten measured values may be outside the limits 27 mm and 33 mm.

```{r, message=FALSE}
cover <- PROB_BASEVAR(
  Id = 2, Name = "cover", DistributionType = "st",
  DistributionParameters = c(30, 2, 3, 10)
)
plan_cover <- PROB_QUALCON(
  Name = "cover_plan", ConformityType = "dim_k||n", ConformityParam = c(10, 1),
  DefectLimitLower = 27, DefectLimitUpper = 33, var = cover
)
plan_cover$getOC()
cover_post <- cover$updateBasevar(plan_cover)
```

The mean value remains unchanged, because the limits are symmetric about the prior mean. The standard deviation decreases.

```{r, fig.width=6, fig.height=4}
cover_post$plot(prior = TRUE)
```

## Limit-state function object (SYS_LSF)

A limit-state function encapsulates the mathematical definition of failure. It can be provided either as an explicit R function or as an arithmetic expression converted internally into a callable closure.

Variables are used as input values for the functions. The naming in the function header (`function(E, R)`) must correspond to its use in the function body (`function(E,R){R-E}`). The naming must also match the name of the variable itself. Therefore, a `list()` with all used `PROB_BASEVAR()` objects is passed to the SYS_LSF() object (note: the objects with their variable names, here var1 and var2, are to be passed, not the variable name itself). In the case of multiple limit state functions being used, the function can be given a name to make it easier to assign in printouts later.
An objective function() can then be defined via the `$func()` property. Alternatively, `$expr()` can be used to assign an `expression(R-E)`.
By means of `$check()` an internal check run is performed, which ensures that all the base variables are complete and suitable.

### Example SYS_LSF

```{r}
var1 <- PROB_BASEVAR(Name = "E", Description = "Effect", DistributionType = "norm", Mean = 10, Sd = 1)
var2 <- PROB_BASEVAR(Name = "R", Description = "Resistance", DistributionType = "norm", Mean = 15, Sd = 1.5)


lsf <- SYS_LSF(vars = list(var1, var2), name = "LSF XY")
lsf$func <- function(E, R) {
  R - E
}
# you can run this check to see if all transformations of the variables worked, all needed data is available and the set of vars fits to the limit state function
lsf$check()
```

The `$check()` method validates that all associated variables are consistent before analysis begins.

## Reliability algorithms (PROB_MACHINE)

Reliability assessment methods are implemented in modular solver objects called `PROB_MACHINE`. Each machine defines one algorithm via its field `$fCall` (e.g., `"FORM"`, `"SORM"`, `"MC_CRUDE"`) together with optional settings stored in `$options`.

* Level 1 (partial safety factors)
* Level 2 (MVFOSM, FORM, SORM)
* Level 3 (MC_CRUDE, MCIS, MC_SubSam)

Levels 2 and 3 can be implemented by means of this package (in order to validate Level 1 procedures). The accuracy of the calculated failure probabilities - and ultimately, reliability indexes - increases through the levels.

A solution-oriented approach is always modelled in a so-called `PROB_MACHINE()`, i.e., an object that stands for a solution algorithm. Over the property `fCall`, the following reliability assessment methods are available:

Level 2 Methods:
* Mean Value First-Order Second Moment (MVFOSM) 
* Advanced First-Order Reliability Method (FORM)
* Second-Order Reliability Method (SORM)

Level 3 Methods:
* Crude Monte-Carlo Simulation (MC_CRUDE)
* Monte-Carlo Simulation with Importance Sampling (MC_IS)
* Monte-Carlo Simulation with Subset Sampling (MC_SubSam)

The property `$name` serves again for better identification in the later output protocol. A list() object in the `$options` property can then be used to specify special options according to the methods.


### MVFOSM

Options:

* h - Accuracy of partial derivatives (finite difference method).
* isExpression - FALSE/TRUE MVFOSM can handle both expressions and functions.

Example:

```{r}
machine <- PROB_MACHINE(name = "MVFOSM", fCall = "MVFOSM", options = list("isExpression" = FALSE, "h" = 0.1))
```
### FORM

Options:

* n_optim - max. amount of evaluation runs
* loctol - local accuracy of the optimisation algorithm (if lower then this threshold, the optim stops)
* optim_type - solver, one of
    * `"auto"` (default) - SLSQP, with improved HLRF (iHLRF) as fallback if SLSQP does not converge. The result field `fallback` shows whether the fallback was used.
    * `"slsqp"` - sequential quadratic programming (package nloptr), no fallback. Good for low-dimensional spaces
    * `"ihlrf"` - improved HLRF with line search, no fallback.
    * `"rackfies"` - Rackwitz-Fiessler iteration (classic HLRF), no fallback.
    * `"auglag"` - Augmented Lagrange solver (nloptr) with random starts.

Example:

```{r}
machine <- PROB_MACHINE(name = "FORM Rack.-Fies.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.001, "optim_type" = "rackfies"))

machine <- PROB_MACHINE(name = "FORM Lagrange.", fCall = "FORM", options = list("optim_type" = "auglag"))
```
### SORM

Options:

* No settings possible

Example:

```{r}
machine <- PROB_MACHINE(name = "SORM", fCall = "SORM")
```
### MC_Crude

Options:

* cov_user=0.05 - abort criterion for the CoV.
* n_batch=10000 - samples per step.
* n_max=1e7 - abort criterion for the maximum number of samples.
* use_threads=NULL - number of workers for backend = "parallel" (default 1). Without effect for backend = "future".
* backend=NULL - "future" or "parallel". If not set, the calculation runs sequentially.
* dataRecord=TRUE - if TRUE, the results per step are returned in `data` (used for visualization of beta).
* debug.level=0 - verbosity level (0 = silent, 1 = summary, 2 = detailed).
* seed - optional integer seed for reproducible random numbers.

Example:

```{r}
machine <- PROB_MACHINE(name = "MC CoV 0.05", fCall = "MC_CRUDE", options = list("n_max" = 1e6, "cov_user" = 0.05))
```

### MC_IS

Options:

* cov_user=0.05 - abort criterion for the CoV.
* n_batch=5000 - samples per step.
* n_max=1e6 - abort criterion for the maximum number of samples.
* use_threads=NULL - number of workers for backend = "parallel" (default 1). Without effect for backend = "future".
* backend=NULL - "future" or "parallel". If not set, the calculation runs sequentially.
* sys_type="serial" - "serial" or "parallel", only used for systems.
* dataRecord=TRUE - if TRUE, the results per step are returned in `data` (single limit-state function only).
* dps - optional design points in physical space. If supplied, FORM is skipped (see `help("MC_IS")` for the format for systems).
* stability_mode="robust" - "robust" or "fast", without influence on the result.
* debug.level=0 - verbosity level (0 = silent, 1 = summary, 2 = detailed).
* seed - optional integer seed for reproducible random numbers.

Further options can be found in the documentation of MC_IS. Changing the options is not recommended for inexperienced users.

```{r, eval = FALSE}
help("MC_IS")
```

Example:

```{r}
machine <- PROB_MACHINE(name = "MC IS", fCall = "MC_IS", options = list("cov_user" = 0.05, "n_max" = 300000))
```
### MC_SubSam

Options:

* Nsubset=1e5 - number of samples in each simulation level.
* p0=0.1 - level probability or conditional probability.
* MaxSubsets=10 - maximum number of simulation levels that are used to terminate the simulation procedure to avoid infinite loop when the target domain cannot be reached.
* Alpha=0.05 - confidence level.
* variance="uniform" - "gaussian" or "uniform".
* debug.level - Integer verbosity level (0 = silent, 1 = summary, 2 = detailed).
* seed - Optional integer seed for reproducible random numbers. 

Example:

```{r}
machine <- PROB_MACHINE(name = "MC Subset Sampling", fCall = "MC_SubSam", options = list("variance" = "uniform"))
```


## System objects

The system object combines all the previous objects and organises the correct calculation of the mathematical model.
Two main classes exist:

* `SYS_PROB` – classical probabilistic system analysis; may contain several LSFs,
* `SYS_PARAM` – parametric sweep over one LSF using varying input parameters.

The `SYS_LSF` object already contains a list of the required basic variables. However, the system object still requires a list of functions to be calculated and the list of algorithms to be used.
In addition, the `$sys_type` property can be used to specify whether the equations have a series `$sys_type="serial"` or parallel `$sys_type="parallel"` relationship. 

The calculation process is then triggered with `$runMachines()`. In the console, the user is informed about the current state of the calculation. After completing the calculation, the results can be retrieved via the `SYS_PROB` object (in the example below `ps`). For more information see result interpretation.

You can set a debug.level field of SYS_PROB to 1 or 2, to get more information during the calculation, by default that value is zero.


### Example – single problem

```{r, results="hide"}
var1 <- PROB_BASEVAR(Name = "E", Description = "Effect", DistributionType = "norm", Mean = 10, Sd = 1)
var2 <- PROB_BASEVAR(Name = "R", Description = "Resistance", DistributionType = "norm", Mean = 15, Sd = 1.5)


lsf <- SYS_LSF(vars = list(var1, var2), name = "LSF XY")
lsf$func <- function(E, R) {
  R - E
}
# you can run this check to see if all transformations of the variables worked, all needed data is available and the set of vars fits to the limit state function
lsf$check()


machine <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.001, "optim_type" = "rackfies"))

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(machine),
  debug.level = 0
)

ps$runMachines()
```

Retrieve results

```{r}
ps$beta_single
```

### Example – system calculations

For systems consisting of multiple LSFs connected either in series or in parallel:

```{r}
# System calculation (System of two identical LSFs)
ps_sys <- SYS_PROB(
  sys_input = list(lsf, lsf),
  probMachines = list(machine),
  sys_type = "serial"
)

ps_sys$runMachines()
ps_sys$calculateSystemProbability()
```

Resulting system reliability indices:

```{r}
ps_sys$beta_sys
```

The method of the system calculation is selected with the argument `calcType` of `$calculateSystemProbability()`:

* `"simpleBounds"` (default) - bounds from the single reliability indices of `$runMachines()`. The machines can be FORM, SORM or other methods.
* `"MCIS"` - importance sampling for serial and parallel systems.
* `"MCC"` - crude Monte Carlo simulation of the system limit-state function.
* `"MCSUS"` - subset simulation of the system limit-state function.

A system analysis by FORM or SORM is not available. These methods only provide the single results for the simple bounds.

## Parametric studies (SYS_PARAM)

Parametric analyses use the class `SYS_PARAM`. This extends `SYS_PROB` to iterate over varying input parameters defined by one or more `PARAM_BASEVAR` objects.

Example:

```{r}
# Force in kN
pvar_F_mod <- PARAM_BASEVAR(
  Name = "Force", DistributionType = "gumbel", Sd = 2, Package = "evd",
  ParamType = "Mean", ParamValues = c(17, 18, 19, 20) 
)

# yield strength in kN, for Mean the mean of the unshifted distribution must be specified.
fy <- PROB_BASEVAR(
  Name = "fy", DistributionType = "slnorm", Package = "TesiproV", Mean = 265000, Sd = 25000, x0 = 160000
)

# plastic moment of resistance in m³ - considered deterministic in this example
W_pl <- 2.5e-4 
# length of cantilever beam in m - considered deterministic in this example
l    <- 2     

lsf_param <- SYS_LSF(vars = list(pvar_F_mod, fy), name = "Spaethe (1992) - Beispiel 2")
lsf_param$func <- function(Force,fy) {
  return(W_pl*fy-l*Force)
}

machine_form <- PROB_MACHINE(name = "FORM", fCall = "FORM")

ps_param <- SYS_PARAM(
  sys_input = list(lsf_param),
  probMachines = list(machine_form)
)

ps_param$runMachines()
print(ps_param$beta_params)
```

Results per iteration can be easily visualised:

```{r}
para_values <- c(17, 18, 19, 20) 
beta_values <- ps_param$beta_params
plot(para_values, beta_values, type = "l", xlab = "Force F [kN]", ylab = "Reliability index β [-]")
```

## Parallelisation in TesiproV
### General
TesiproV supports parallel execution for the Monte Carlo-based reliability methods `MC_IS()` and `MC_CRUDE()`. `MC_SubSam()` runs in a single process. The package provides two parallel backends:

* "parallel" -- based on base R's parallel package  
* "future" -- based on the future ecosystem

`MC_IS()` does not modify the active future plan, users configure parallel execution explicitly. `MC_IS()` and `MC_CRUDE()` do not modify the active future plan, users configure parallel execution explicitly.

### Reproducibility and Determinism
`MC_IS()` is designed to be:

* Deterministic for a fixed seed
* Independent of the number of worker threads
* Independent of the selected backend

This means that:
```{r, eval=FALSE}
MC_IS(..., seed = 123, use_threads = 1)
MC_IS(..., seed = 123, use_threads = 8)
```

will produce identical results up to rounding errors.

Determinism is guaranteed by performing all random number generation in the master process. Worker processes do not generate random numbers.

`MC_CRUDE()` follows the same principle: all random numbers are generated in the master process, so the result for a fixed seed does not depend on the number of workers.

### Backend Overview
1. The "parallel" Backend

* Uses `parallel::mclapply()` on Linux/macOS  
* Uses PSOCK clusters on Windows  
* Suitable for high-performance Linux clusters  
* No dependency on `future`

Example:

```{r, eval=FALSE}
MC_IS(..., backend = "parallel", use_threads = 8)
```

2. The "future" Backend
   
Delegates parallelisation to the future framework.
TesiproV does not change the global future plan. Users must configure parallel execution explicitly.

Linux / macOS example:
```{r, eval=FALSE}
library(future)
plan(multicore, workers = 8)

MC_IS(..., backend = "future")
```

Windows example:
```{r, eval=FALSE}
library(future)
plan(multisession, workers = 8)

MC_IS(..., backend = "future")
# If no future plan is set, execution will be sequential.
```

### Recommended Parallel Settings

Parallel efficiency depends strongly on the computational cost of the limit-state function (LSF).

| Problem Type               | Recommended Threads | Backend Recommendation      | Suggested `n_batch`      |
|----------------------------|--------------------|-----------------------------|--------------------------|
| Very simple LSF           | 1-2                | `"future"` or `"parallel"`  | 50,000 - 100,000         |
| Moderately expensive LSF  | 4-8                | `future(multicore)`         | 100,000 - 250,000        |
| Expensive LSF             | 8-16               | `future(multicore)`         | 250,000 - 500,000        |
| Very expensive LSF        | ~16                | `future(multicore)`         | ~500,000                 |


Using more than 16 threads on large shared-memory systems often leads to memory-bandwidth saturation.

### Batch Size Selection
Performance improves significantly when n_batch is large enough to amortize parallel overhead.

Typical observations:

* Small batches -> overhead dominates
* Medium batches -> optimal performance
* Very large batches -> memory-bound regime

For HPC workloads, a batch size between:

```{r, eval=FALSE}
250,000 -- 500,000
```

often provides near-optimal performance.

### Performance Characteristics
Parallel scaling depends on:

* LSF computational complexity
* Memory bandwidth
* NUMA architecture
* Batch size

For CPU-intensive LSFs, typical scaling on modern multi-core systems is:

* 2 threads -> ~1.6x speedup  
* 4 threads -> ~2.8-3.5x  
* 8 threads -> ~4-6×  
* 16 threads -> ~6-9×  
* 32 threads -> often limited by memory bandwidth


### HPC Usage Notes
On large Linux clusters (e.g., NUMA systems):

* `future(multicore)` typically provides the best performance  
* Avoid very small n_batch values  
* Prefer 8-16 threads for optimal efficiency  
* Use deterministic seeds for reproducible studies

Example HPC workflow:
```{r, eval=FALSE}
library(future)
plan(multicore, workers = 16)

res <- MC_IS(
  lsf = my_lsf,
  lDistr = my_distr,
  n_batch = 300000,
  backend = "future",
  seed = 12345
)
```

## Interpretation of results and post processing

After computation each system object contains structured result lists:
For regular problems (`SYS_PROB`):

* Single calculations -> `$beta_single`, `$res_single`,
* System calculations -> `$beta_sys`, `$res_sys`.

Here beta_* is the beta values summarised while res_* is the full result output via list(). The list `$res_*` is structured as follows `ps$res_*[[lsf]][[machine]]$value`.

For parametric studies (`SYS_PARAM`):

* Parameter iterations -> `$beta_params`, `$res_params`.

where the first list level of `$res_params` corresponds to the parameter iterations. The data can be then processed individually.

In addition, the results can be printed to a file using the command `$printResults(filename)`, where filename corresponds to the file name and the file is stored in the current working directory (see `getwd()` or `setwd()`).

The project can also be stored differently to allow later postprocessing in other files. For this, the command `$saveProject(level,filename)` is used. The following stages describe a different level of detail:

* 1. Only the beta values are stored.
* 2. The `$res_*` object is saved.
* 3. The entire `ps` object is stored.
* 4. An image of the current working memory is created, which can be reloaded later.

Stages 1 to 3 can be then read back into working memory using the `readRDS(filename)` command, while Stage 4 is reloaded using `load(filename)`.

Example:

```{r, eval = FALSE}
ps <- SYS_PARAM(
  sys_input = list(lsf),
  probMachines = list(machine)
)

ps$runMachines()

local_beta <- ps$beta_params
local_runtime <- ps$res_params[[1]][[1]][[1]]$runtime

# define path with setwd("...")
ps$printResults("result_file") # report file result_file.txt
ps$saveProject(2, "project_file") # results in project_file_res.rds
ps$saveProject(4, "project_file") # workspace in project_file.RData
```

The results can be reloaded later via:

```{r, eval = FALSE}
readRDS("project_file_res.rds")
load("project_file.RData")
```

Debug output verbosity is controlled by field $debug.level (0=silent … 2=detailed).

# Examples

## Two normally distributed variables
A limit state equation 20-(x1-x2)^2-8*(x1+x2-4)^3 is to be investigated, where X1 and X2 are normally distributed around the expected value of 0.25 and the standard deviation of 1. The example is taken from the documentation of the UQLab software. 

Firstly, two basic variables are created and assigned to the R variables Var1 and Var2. The first variable is named as "X1" and the second as "X2". These names are then reused in the definition of the limit state function (i.e., to compare the input variables of the function). In addition, both variables receive the indicator "norm" as distribution type (DistributionType) and the mean value defined above or the corresponding standard deviation.

Subsequently, a limit state function object (SYS_LSF) is assigned to the R variable lsf. A list of variables (Var1,Var2) is passed to the object and a name is assigned for later identification. In addition, the corresponding limit state equation is defined via the object property "func".

The R variable "FORM" is assigned a solution method of the Level 2 (FORM) without any additional parameters.

Finally, the modules are merged in the SYS_PROB object and the R variable "ps" and the calculation is started by means of the $runMachines() method.
In the following result log the intermediate results of the individual iteration steps are shown. Using the "ps" object, different results (here, for example, the beta value) can be retrieved after a successful calculation.


```{r, results="hide"}
Var1 <- PROB_BASEVAR(Id = 1, Name = "X1", DistributionType = "norm", Mean = 0.25, Sd = 1) # kN/m²
Var2 <- PROB_BASEVAR(Id = 2, Name = "X2", DistributionType = "norm", Mean = 0.25, Sd = 1) # m

lsf <- SYS_LSF(vars = list(Var1, Var2), name = "UQLab 2d_hat")
lsf$func <- function(X1, X2) {
  return(20 - (X1 - X2)^2 - 8 * (X1 + X2 - 4)^3)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(form)
)
ps$runMachines()
```

```{r}
ps$beta_single
```

## Non-normally distributed basic variable
To demonstrate non-normally distributed basic variables, an example from Nowak & Collins (Example 5.11, Page 127/128) is calculated. Here normally distributed, lognormally distributed as well as gumbel distributed variables are used. While "norm" and "lnorm" are still contained in the standard scope of R (package "stats" belongs to base R) the Gumbel distribution is to be considered on the basis of the EVD-package (EVD := extreme value distributions).

```{r, results ='hide'}
X1 <- PROB_BASEVAR(Id = 1, Name = "Z", DistributionType = "norm", Mean = 100, Cov = 0.04) # kN/m²
X2 <- PROB_BASEVAR(Id = 2, Name = "Fy", DistributionType = "lnorm", Mean = 40, Cov = 0.1) # m
X3 <- PROB_BASEVAR(Id = 2, Name = "M", DistributionType = "gumbel", Package = "evd", Mean = 2000, Cov = 0.1) # m

lsf <- SYS_LSF(vars = list(X1, X2, X3), name = "Nowak & Collins Exp5.11_1")
lsf$func <- function(M, Fy, Z) {
  return(Z * Fy - M)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf),
  probMachines = list(form)
)
ps$runMachines()
```

```{r}
ps$beta_single
```



## Multiple limit state functions in one system

Basically, the "SYS_PROB" object can calculate multiple LSF's using both Level 2 and Level 3 methods. This can be used when a system consists of several limit state functions and when the system reliability is to be determined using SimpleBounds. In this case, the failure probability must first be determined for each individual limit state equation. These results are then used to calculate bounds for the system failure probability. On the other hand, the development of a limit state function often poses several variants in equations or variables, so that by the simultaneous calculation of several equations satisfactory comparability is achieved and a high usability remains granted.

However, if the system failure probability can not to be sufficiently estimated by means of the SimpleBounds, simulation methods (e.g., MC-Crude) can be used to calculate the system failure probability more exactly. Then, the system can also be constructed according to the following scheme.

In this example, the above examples are first determined individually as a system using a FORM calculation and then the system failure probability is determined using different methods.



```{r, results="hide"}
Var1 <- PROB_BASEVAR(Id = 1, Name = "X1", DistributionType = "norm", Mean = 0.25, Sd = 1) # kN/m²
Var2 <- PROB_BASEVAR(Id = 2, Name = "X2", DistributionType = "norm", Mean = 0.25, Sd = 1) # m

lsf1 <- SYS_LSF(vars = list(Var1, Var2), name = "UQLab 2d_hat")
lsf1$func <- function(X1, X2) {
  return(20 - (X1 - X2)^2 - 8 * (X1 + X2 - 4)^3)
}


X1 <- PROB_BASEVAR(Id = 4, Name = "Z", DistributionType = "norm", Mean = 100, Cov = 0.04) # kN/m²
X2 <- PROB_BASEVAR(Id = 5, Name = "Fy", DistributionType = "lnorm", Mean = 40, Cov = 0.1) # m
X3 <- PROB_BASEVAR(Id = 6, Name = "M", DistributionType = "gumbel", Package = "evd", Mean = 2000, Cov = 0.1) # m

lsf2 <- SYS_LSF(vars = list(X1, X2, X3), name = "Nowak & Collins Exp5.11_1")
lsf2$func <- function(M, Fy, Z) {
  return(Z * Fy - M)
}

form <- PROB_MACHINE(name = "FORM Rack.-Fieß.", fCall = "FORM")

ps <- SYS_PROB(
  sys_input = list(lsf1, lsf2),
  sys_type = "serial",
  probMachines = list(form)
)
ps$runMachines()
```

```{r}
ps$beta_single

ps$calculateSystemProbability()
ps$beta_sys
```


The following lines are not executed because the calculation time is several minutes depending on the machine.
Results in the range of the above interval are to be expected. The results can be retrieved after the calculation using "$beta_sys".
```{r, eval = FALSE}
params <- list("cov_user" = 0.05, "n_max" = 50000)
ps$calculateSystemProbability("MCC", params)

ps$calculateSystemProbability("MCIS", params)

params_sus <- list("Nsubset" = 10000, "p0" = 0.1)
ps$calculateSystemProbability("MCSUS", params_sus)
```

## Reliability with and without quality control

The posterior variable can be used like any other basic variable. In this example the concrete strength `fc` from the section on quality control is the resistance and a normally distributed action `E` (mean 25, standard deviation 3, same unit) is the load effect. The reliability index is calculated with the prior and with the posterior distribution of the strength.

```{r, results="hide", message=FALSE}
E <- PROB_BASEVAR(Id = 3, Name = "E", Description = "Action effect",
                  DistributionType = "norm", Mean = 25, Sd = 3)

form <- PROB_MACHINE(name = "FORM", fCall = "FORM")

lsf_prior <- SYS_LSF(vars = list(fc, E), name = "Prior")
lsf_prior$func <- function(fc, E) fc - E
lsf_prior$check()

lsf_post <- SYS_LSF(vars = list(fc_post, E), name = "Posterior")
lsf_post$func <- function(fc, E) fc - E
lsf_post$check()

ps <- SYS_PROB(sys_input = list(lsf_prior, lsf_post), probMachines = list(form))
ps$runMachines()
```

```{r}
ps$beta_single
```

The reliability index is higher after the acceptance of the lot. This is the effect of the quality control on the safety level.

## Parameterstudy

* Punching according to DIN EN 1992-1 + NA(D)
* Effective depth (var.d) is parameterised (cf. vector at "ParamValues")
* Additional parametric deterministic variable "var.V_Ed" which depends on "pre.d
* Calculation by FORM and postprocessing see below.


```{r, results="hide", warning=FALSE}
var.d <- PARAM_BASEVAR(Id = 1, Name = "d", Description = "Effective depth", DistributionType = "norm", Sd = 10, ParamType = "Mean", ParamValues = c(seq(200, 800, 50), seq(1000, 3000, 100)) + 10)
pre.d <- c(seq(200, 800, 50), seq(1000, 3000, 100))

var.f_ck <- PROB_BASEVAR(Id = 2, Name = "f_ck", Description = "Concrete compressive strength", Package = "TesiproV", DistributionType = "lt", DistributionParameters = c(3.85, 0.09, 3, 10))
pre.f_ck <- 35

var.f_ywk <- PROB_BASEVAR(Id = 3, Name = "f_ywk", Description = "Yield strength of shear reinforcement", DistributionType = "norm", Mean = 560, Cov = 0.02)
pre.f_ywk <- 500

var.f_yk <- PROB_BASEVAR(Id = 4, Name = "f_yk", Description = "Yield strength of longitudinal reinforcement", DistributionType = "norm", Mean = 560, Sd = 30)
pre.f_yk <- 500

var.rho_l <- PROB_BASEVAR(Id = 5, Name = "rho_l", Description = "Longitudinal reinforcement ratio", DistributionType = "norm", Mean = 0.01, Cov = 0.02)
pre.rho_l <- var.rho_l$Mean

var.c_D <- PROB_BASEVAR(Id = 6, Name = "c_D", Description = "Column diameter", DistributionType = "norm", Mean = 500 + 0.003 * 500, Sd = 4 + 0.006 * 500)
pre.c_D <- 500

var.theta_c <- PROB_BASEVAR(Id = 7, Name = "theta_c", Description = "Model uncertainty", DistributionType = "lnorm", Mean = 1.0644, Cov = 0.1679)

var.gamma_c <- PROB_DETVAR(Id = 8, Name = "gamma_c", Description = "Partial safety factor concrete", Value = 1.5)
var.gamma_s <- PROB_DETVAR(Id = 9, Name = "gamma_s", Description = "Partial safety factor steel", Value = 1.15)
var.alpha_cc <- PROB_DETVAR(Id = 10, Name = "alpha_cc", Description = "Long-term factor alpha_cc", Value = 0.85)

V_Ed <- c(
  0.6412447, 0.8760515, 1.1424185, 1.4399554, 1.7554079, 2.0188587, 2.3003489, 2.5997038, 2.9167825, 3.2514671, 3.6036567, 3.9732630,
  4.3602076, 6.0800504, 7.0424326, 8.0726216, 9.1702994, 10.3351829, 11.5670173, 12.8655715, 14.2306347, 15.6620132, 17.1595283, 18.7230145,
  20.3523175, 22.0472935, 23.8078076, 25.6337332, 27.5249511, 29.4813487, 31.5028193, 33.5892621, 35.7405811, 37.9566849
)
var.V_Ed <- PARAM_DETVAR(Id = 11, Name = "V_Ed", Description = "Applied shear force", ParamValues = V_Ed)

lsf1 <- SYS_LSF(
  vars = list(
    var.d, var.f_ck, var.f_yk, var.f_ywk, var.rho_l, var.c_D,
    var.theta_c, var.V_Ed, var.gamma_c, var.gamma_s, var.alpha_cc
  ),
  name = "Durchstanzwiderstand ohne Bewehrung"
)
lsf1$func <- function(d, f_ck, f_yk, f_ywk, rho_l, c_D, theta_c, V_Ed, gamma_c, gamma_s, alpha_cc) {
  f_cd <- f_ck * alpha_cc / gamma_c
  f_yd <- f_yk / gamma_s
  f_ywd <- f_ywk / gamma_s

  u_0 <- c_D * pi
  u_1 <- (c_D / 2 + 2 * d) * 2 * pi
  u_0.5 <- (c_D / 2 + 0.5 * d) * 2 * pi

  k <- min(1 + sqrt(200 / d), 2)

  if ((u_0 / d) < 4) {
    C_Rd_c <- 0.18 * (0.1 * u_0 / d + 0.6)
  } else {
    C_Rd_c <- 0.18
  }
  rho_l <- min(rho_l, 0.5 * f_cd / f_yd, 0.02)
  C_Rd_c1 <- C_Rd_c * k * (100 * rho_l * f_ck)^(1 / 3)

  v_Rd_c <- C_Rd_c1
  V_Rd_c <- v_Rd_c * d * u_1 * 1e-6 * theta_c
  return(V_Rd_c - V_Ed)
}

form <- PROB_MACHINE(name = "FORM RF.", fCall = "FORM", options = list("n_optim" = 20, "loctol" = 0.0001, "optim_type" = "rackfies"))

ps <- SYS_PARAM(
  sys_input = list(lsf1),
  probMachines = list(form)
)
ps$runMachines()
```

```{r, results="hide", warning=FALSE}
d <- pre.d
beta <- unlist(unname(ps$beta_params))
plot(x = d, y = beta, type = "l")
```


# References

1. JOINT COMMITTEE ON STRUCTURAL SAFETY, JCSS Probabilistic Model Code. URL: https://www.jcss-lc.org/, 2001. 

2. Marelli & Sudret (2014): UQLab: A Framework for Uncertainty Quantification in MATLAB, ICVRAM 2014 Conference Proceedings. DOI:10.1061/9780784413609.257

3. NOWAK, A. S.; COLLINS, K. R. Reliability of structures. CRC Press, 2012.

4. DIN EN 1992-1-1: 2011-01: DIN EN 1992-1-1: Design of concrete structures – Part 1-1: General rules and rules for buildings; German version EN 1992-1-1:2004 + AC:2010. 

5. DIN EN 1992-1-1/NA: 2013-04: National Annex – Nationally determined parameters – DIN EN 1992-1-1: Design of concrete structures – Part 1-1: General rules and rules for buildings. 

# Acknowledgement

The TesiproV package was developed within the framework of the research project “TesiproV”, funded by the German Federal Ministry for Economic Affairs and Energy (BMWi), and the research project “Qualitätskontrollen und Bauwerkszuverlässigkeit im Bauwesen”, funded by the German Institute for Building Technology (DIBt). The authors gratefully acknowledge the financial support provided by these institutions.

# Contact & License

For questions, feedback or collaboration requests please contact us at:

**TU Dortmund University – Chair of Structural Concrete**  
August-Schmidt-Strasse 8 · D‑44227 Dortmund · Germany  

Email: [tesiprov.lmb.ab@tu-dortmund.de](mailto:tesiprov.lmb.ab@tu-dortmund.de)  

Additional contacts:  
* [nille-hauf@hochschule-bc.de](mailto:nille-hauf@hochschule-bc.de) ([konstantin.nille-hauf@googlemail.com](mailto:konstantin.nille-hauf@googlemail.com))  
* [tania.feiri@tu-dortmund.de](mailto:tania.feiri@tu-dortmund.de)  
* [marcus.ricker@tu-dortmund.de](mailto:marcus.ricker@tu-dortmund.de)  
* [til.lux@tu-dortmund.de](mailto:til.lux@tu-dortmund.de)

---

© 2021–2026 K. Nille-Hauf, T. Feiri, M. Ricker, T. Lux -- Hochschule Biberach (until 2022), TU Dortmund University – Chair of Structural Concrete (since 2023)

TesiproV is released under the [MIT License](https://opensource.org/licenses/MIT).  
You are free to use, modify and redistribute the software provided that this copyright notice and license information remain included.

For official releases see the [CRAN package TesiproV](https://CRAN.R-project.org/package=TesiproV).  
Source code repository: <https://gitlab.tu-dortmund.de/ls-massivbau/tesiprov>