#' @title First-order reliability method
#'
#' @description The First-Order Reliability Method computes an estimation of the failure
#' probability by approximating the limit-state function at the Most Probable Failure Point
#' with a hyperplane.
#'
#' @details The FORM method has to be used in the standard Gaussian input space. It is designed
#' to estimate probability of the form \eqn{P[g(\mathbf{X}) < 0]}{P[g(X) < 0]} with g the limit-state function.
#' This function has to be modified accordingly to fit into this framework
#'
#' Furthermore, it should be able to handle matrix input of column vectors. See the mistral vignette
#' for more info about \code{lsf} definition
#'
#' @author Vincent MOUTOUSSAMY and Clement WALTER \email{clementwalter@icloud.com}
#'
#' @references
#' \itemize{
#' \item
#' O. Ditlevsen and H.O. Madsen. Structural reliability methods, Wiley, 1996\cr
#' \item
#' M. Lemaire, A. Chateauneuf and J. Mitteau. Structural reliability, Wiley Online Library, 2009.\cr
#' }
#'
#' @examples
#' \dontrun{
#' # u.dep is a starting point for the research of the Most Probable Failing Point
#' # N.calls is a total number of calls
#' form <- mistral::FORM(dimension = 2, mistral::kiureghian, N.calls = 1000,
#' u.dep = c(0,0))
#' form$p
#'
#' # use IS=TRUE to use an Importance Sampling scheme with a Gaussian standard
#' # proposal distribution centred at the MPFP
#' form.IS <- mistral::FORM(dimension = 2, mistral::kiureghian, N.calls = 1000,
#' u.dep = c(0,0),
#' IS = TRUE)
#' form.IS$p
#' }
#'
#' @return A list containing the following objects
#' \item{p}{Failure probability}
#' \item{indice.reliab}{Reliability index}
#' \item{Ncall}{Number of calls to f}
#' \item{Design.Point}{Coordinates of the design point}
#' \item{fact.imp}{Importance factors}
#' \item{variance}{Standard error of the probability estimator (if IS = TRUE)}
#' \item{Interval.conf}{Confidence interval of the estimator at 0.95 (if IS = TRUE)}
#' \item{DOE}{List which contains the design of experiments}
#'
#' @import ggplot2
#'
#' @export
FORM <- function(dimension,
#' @param dimension the dimension of the input space.
lsf,
#' @param lsf the limit-state function.
u.dep = rep(0, dimension),
#' @param u.dep the starting point for the MPFP search.
N.calls = 1e2,
#' @param N.calls the total number of calls for the whole algorithm.
eps = 1e-7,
#' @param eps stopping criterion: distance of two points between two iterations.
Method = "HLRF",
#' @param Method choice of the method to search the design point: "AR" for Abdo-Rackwitz
#' and "HLRF" for Hasofer-Lindt-Rackwitz-Fiessler.
IS = FALSE,
#' @param IS "TRUE" for using importance Sampling method with an standard Gaussian importance
#' density centred at the MPFP.
IS.ratio = 0.5,
#' @param IS.ratio ratio of N.calls for the search of the design point by FORM. Default = 0.5.
#' 1-\code{IS.ratio} = the remaining ratio to be used for importance sampling.
plot = FALSE,
#' @param plot to plot the generated samples.
plot.lsf = FALSE,
#' @param plot.lsf a boolean indicating if the \code{lsf} should be added to the
#' plot. This requires the evaluation of the \code{lsf} over a grid and
#' consequently should be used only for illustation purposes.
plot.lab = c('x_1', 'x_2')
#' @param plot.lab the x and y labels for the plot.
){
# Few error messages :
if (mode(lsf) != "function") stop("lsf is not a function")
if(plot==TRUE & dimension>2){
message("Cannot plot dimension>2 => plot <- FALSE")
plot <- FALSE
}
# Fix check NOTE
x <- y <- z <- NULL
#------ Gradient------#
GRAD <- function(func, X, f.X){
X <- as.matrix(X)
dX <- X[,rep(1, dimension)] + 1e-6*diag(dimension)
f.dX <- lsf(dX)
dfunc <- (f.dX - f.X)/1e-6
return(list(g = dfunc,
dX = dX,
f.dX = f.dX))
}
#----------------------------------------------------------------------#
List.res <- list() # List of the output
u.dep <- as.matrix(u.dep)
dimnames(u.dep) = list(rep(c('x', 'y'), length.out = dimension))
u.new <- 0
cp <- 0 #Numbers of calls to f
fact.imp <- 0 #importance factor
N.FORM <- ifelse(IS, floor(N.calls*IS.ratio), N.calls) #Number of calls reserved to FORM
u.ans <- lsf(u.dep); cp <- cp + 1
DOE <- list(x = u.dep, y = u.ans) # Design of experiments
if(plot==TRUE){
xplot <- yplot <- c(-80:80)/10
df_plot = data.frame(expand.grid(x=xplot, y=yplot))
df.points <- data.frame(t(DOE$x), z = DOE$y)
p <- ggplot(data = df.points, aes(x,y)) +
theme(legend.position = "none") +
xlim(-8, 8) + ylim(-8, 8) + xlab(plot.lab[1]) + ylab(plot.lab[2]) +
geom_point(aes(color=z))
if(plot.lsf==TRUE){
zplot = lsf(t(df_plot))
p <- p + geom_contour(data = data.frame(df_plot, z = zplot), aes(z=z, color=z), breaks = 0)
}
print(p)
}
res <- GRAD(lsf, u.dep, u.ans); cp <- cp + dimension
g <- res[[1]]
DOE$x <- cbind(DOE$x, res$dX)
DOE$y <- c(DOE$y, res$f.dX)
err <- 1 #distance of two tested points
#----------------------------------Algorithme HLRF :----------------- #
if(Method == "HLRF"){
a <- g/sqrt(sum(g^2))
while( (err > eps) & (cp <= (N.FORM - (1 + 2*dimension)) ) ){ # modif BIS
u.new <- (sum(u.dep*a)*a - u.ans*a/sqrt(sum(g^2)))
err <- sqrt(sum((u.dep-u.new)^2))
u.dep <- u.new
u.ans <- lsf(u.dep); cp <- cp + 1
DOE$x <- cbind(DOE$ x, u.dep)
DOE$y <- c(DOE$y, u.ans)
res <- GRAD(lsf, u.dep, u.ans); cp <- cp + dimension
g <- res[[1]]
DOE$x <- cbind(DOE$x, res$dX)
DOE$y <- c(DOE$y, res$f.dX)
a <- g*ifelse( sum(g^2) == 0, 0, 1/sqrt(sum(g^2)))
}
}
#-----------------------------End of the algorithm HLRF :--------------------------#
#-----------------------------Algorithm Abdo-Rackwitz :--------------------------#
if(Method == "AR"){
lambda <- 2*(u.ans - (t(g)%*%u.dep))/(g%*%g)
while( (err > eps) & (cp <= (N.FORM - (1 + 2*dimension)) ) ){ # modif BIS
u.new <- -lambda*g/2
err <- sqrt(t(u.dep-u.new)%*%(u.dep-u.new))
u.dep <- u.new
u.ans <- lsf(u.dep); cp <- cp + 1
DOE$x <- cbind(DOE$ x, u.dep)
DOE$y <- c(DOE$y, u.ans)
res <- GRAD(lsf, u.dep, u.ans); cp <- cp + dimension
g <- res[[1]]
DOE$x <- cbind(DOE$x, res$dX)
DOE$y <- c(DOE$y, res$f.dX)
lambda <- 2*(u.ans - sum(g*u.dep))/sum(g^2)
if(cp >= N.FORM){break}
}
}
#------------------------End of the algorithm of Abdo-Rackwitz----------------------#
if(plot==TRUE){
p$data <- data.frame(t(DOE$x), z = DOE$y, row.names = NULL, check.names = FALSE)
print(p)
}
a <- g/sqrt(sum(g^2))
B <- -a%*%u.dep
P <- c(pnorm(-B))
fact.imp <- a^2
List.res <- list(P, B, cp, u.dep, fact.imp, DOE);
names(List.res) <- c("p", "indice.reliab", "Ncall", "Design.Point", "fact.imp", "DOE")
#---------------------------Importance Sampling :---------------------------------#
if(IS == TRUE){
alpha <- 0.05
N.IS <- N.calls - cp #Number of call reserved to the importance sampling
norm.u <- u.dep%*%u.dep
P.temp <- 0
P <- 0
VAR <- 0
s <- 0
Z <- numeric(N.IS)
v.temp <- matrix(rnorm(N.IS*dimension,0,1),ncol=dimension)
v.temp <- v.temp + matrix(rep(u.dep,N.IS),ncol=dimension,byrow=TRUE) #One simulate points around u.dep
zz <- v.temp - matrix(rep(u.dep,N.IS), ncol = dimension, byrow = TRUE)
s <- exp( u.dep%*%u.dep/2 - as.numeric(v.temp%*%u.dep) )
test.v.temp <- lsf(t(v.temp))
DOE$x <- cbind(DOE$x, t(v.temp))
DOE$y <- c(DOE$y, test.v.temp)
indic.Def <- 1*(test.v.temp <=0)
P.temp <- indic.Def*s
List.res$Ncall <- cp + N.IS #Numbers of call to the failure fonction
List.res$p <- P <- c(sum(P.temp)/N.IS) #failure probability
List.res$indice.reliab <- -qnorm(P) #indices of reliability
List.res$variance <- VAR <- (sum(P.temp^2)/N.IS - P^2)/(N.IS-1) #standard error of the estimator
IC.inf <- P - qnorm(1-alpha/2)*sqrt(VAR) #low bound of confidence interval
IC.sup <- P + qnorm(1-alpha/2)*sqrt(VAR) #high bound of confidence interval
List.res$Interval.conf <- c(IC.inf,IC.sup)
List.res$fact.imp <- a^2
List.res$DOE <- DOE
}
#----------------------------End Importance Sampling----------------------------#
return(List.res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.