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.

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:

Additional placeholders exist for future extensions.

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

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).

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.

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

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.

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.

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()

plan$plot()

fc_post <- fc$updateBasevar(plan)
#> updateBasevar complete.
#>   Prior     : Mean = 42.9141  Sd = 5.8751  (numerical approximation)
#>   Posterior : Mean = 43.3318  Sd = 5.4219

c(Mean = fc_post$Mean, Sd = fc_post$Sd)
#>      Mean        Sd 
#> 43.331813  5.421924

fc_post$plot(prior = TRUE)

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)

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.

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()

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.

machine <- PROB_MACHINE(name = "MVFOSM", fCall = "MVFOSM", options = list("isExpression" = FALSE, "h" = 0.1))

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"))

machine <- PROB_MACHINE(name = "SORM", fCall = "SORM")

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

help("MC_IS")

machine <- PROB_MACHINE(name = "MC IS", fCall = "MC_IS", options = list("cov_user" = 0.05, "n_max" = 300000))

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.

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()

ps$beta_single
#>                    LSF XY
#> FORM Rack.-Fieß. 2.773501

# 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()

ps_sys$beta_sys
#>                               [,1]
#> SB min (FORM Rack.-Fieß.) 2.540293
#> SB max (FORM Rack.-Fieß.) 2.773501

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:

# 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)
#> $`Paramrun No. 1. Paramname Force with value(Type): 17.000000 (Mean)`
#>      Spaethe (1992) - Beispiel 2
#> FORM                    3.775581
#> 
#> $`Paramrun No. 2. Paramname Force with value(Type): 18.000000 (Mean)`
#>      Spaethe (1992) - Beispiel 2
#> FORM                     3.61394
#> 
#> $`Paramrun No. 3. Paramname Force with value(Type): 19.000000 (Mean)`
#>      Spaethe (1992) - Beispiel 2
#> FORM                    3.445885
#> 
#> $`Paramrun No. 4. Paramname Force with value(Type): 20.000000 (Mean)`
#>      Spaethe (1992) - Beispiel 2
#> FORM                    3.270672

Results per iteration can be easily visualised:

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

MC_IS(..., seed = 123, use_threads = 1)
MC_IS(..., seed = 123, use_threads = 8)

MC_IS(..., backend = "parallel", use_threads = 8)

library(future)
plan(multicore, workers = 8)

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

library(future)
plan(multisession, workers = 8)

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

250,000 -- 500,000

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:

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:

readRDS("project_file_res.rds")
load("project_file.RData")

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

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.

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()

ps$beta_single
#>                  UQLab 2d_hat
#> FORM Rack.-Fieß.     3.434565

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).

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()

ps$beta_single
#>                  Nowak & Collins Exp5.11_1
#> FORM Rack.-Fieß.                  4.022115

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.

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()

ps$beta_single
#>                  UQLab 2d_hat Nowak & Collins Exp5.11_1
#> FORM Rack.-Fieß.     3.434565                  4.022115

ps$calculateSystemProbability()
ps$beta_sys
#>                               [,1]
#> SB min (FORM Rack.-Fieß.) 3.409355
#> SB max (FORM Rack.-Fieß.) 3.434565

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”.

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.

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()

ps$beta_single
#>         Prior Posterior
#> FORM 2.893513  3.155135

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.

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()

d <- pre.d
beta <- unlist(unname(ps$beta_params))
plot(x = d, y = beta, type = "l")