2 Tank Rel(Fréchet)

Post date: Apr 01, 2014 9:52:17 AM

#################################################################################

# Compute reliability of Vesely's pressurized tank system from sparse sample data

# without making any assumptions about dependence among the components

# wipe everything from memory

rm(list=ls())

many = 10000

precise <- function(b) if (length(b)==many) TRUE else FALSE

leftside <- function(b) if (precise(b)) return(b) else return(b[1:many])

rightside <- function(b) if (precise(b)) return(b) else return(b[(many+1):(2*many)])

sampleindices = function() round(runif(many)*(many-1) + 1)

env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error') # only works for precise distributional inputs

beta <- function(v,w) if ((v==0) && (w==0)) env(rep(0,many),rep(1,many)) else if (v==0) rep(0,many) else if (w==0) rep(1,many) else sort(rbeta(many, v, w))

pairsides <- function(b) {i = sampleindices(); return(env(leftside(b)[i],rightside(b)[i]))}

# c-box for Bernoulli parameter:  x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)

kn <- function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k)))) 

# the Frechet logical operators, where x = P(X), y = P(Y)

# P(X & Y) = [max(0, x + y – 1), min(x, y)]

# P(X or Y) = [max(x, y), min(1, x + y)]

andF <- function(x,y) env(pmax(0,leftside(x)+leftside(y)-1), pmin(rightside(x),rightside(y)))

orF <- function(x,y) env(pmax(leftside(x), leftside(y)), pmin(1, rightside(x)+rightside(y)))

plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Probability',...) {

  edf <- function (x, col, lwd, ...) {

      n <- length(x)

      s <- sort(x)

      if (2000<n) {s = c(s[1],s[1:(n/5)*5],s[n]); n <- length(s);} # subsetting for the plot to every 5th value

      lines(c(s[[1]],s[[1]]),c(0,1/n),lwd=lwd,col=col,...)

      for (i in 2:n) lines(c(s[[i-1]],s[[i]],s[[i]]),c(i-1,i-1,i)/n,col=col,lwd=lwd,...)

      }

  b = ifelse(b==-Inf, xlim[1] - 10, b)

  b = ifelse(b==Inf, xlim[2] + 10, b)

  if (new) plot(NULL, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

  if (length(b) < many) edf(b,col,lwd) else

  edf(c(min(b),max(b),b[1:min(length(b),many)]),col,lwd)

  if (many < length(b)) edf(c(min(b),max(b),b[(many+1):length(b)]),col,lwd)

  }

#################################################################################

#################################################################################

#################################################################################

# example of c-boxes with Vesely's risk of the tank rupturing under pressurization

################

# hypothetical data in the form of k failures observed out of n trials are used to create c-boxes using the formula kn(k,n)

t = kn(0, 2000)

k2 = kn(3, 500)

s = kn(1, 150)

s1 = kn(0, 460)

k1 = kn(7, 10000)

r = kn(0, 380)

################

# compute failure probability of top event (tank rupturing under pressurization)

# without making any assumptions about the dependence among components

#e2 = t %|% (k2 %|% (s %&% (s1 %|% (k1 %|% r))))

e1 = orF(t, orF(k2, andF(s, orF(s1, orF(k1, r)))))

par(mfrow=c(2,4))

plotbox(t); title('tank T')

plotbox(k2); title('relay K2')

plotbox(s); title('pressure switch S')

plotbox(s1); title('on-off switch S1')

plotbox(k1); title('relay K1')

plotbox(r); title('timer relay R')

plotbox(e1); title('top event E1')

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// check the calculation in Risk Calc

func Beta() if (($1==0) & ($2==0)) return [0,1] else if ($1==0) return 0 else if ($2==0) return 1 else return beta($1,$2)

// c-box for Bernoulli parameter:  x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)

function kn() return((env(Beta($1, $2-$1+1), Beta($1+1, $2-$1))))

t = kn(0, 2000)

k2 = kn(3, 500)

s = kn(1, 150)

s1 = kn(0, 460)

k1 = kn(7, 10000)

r = kn(0, 380)

# logical operators

orI <- function(x,y) return(1-(1-x)*(1-y))

andI <- function(x,y) return(x*y)

# assuming mutual independence

e1 = t ||| (k2 ||| (s |&| (s1 ||| (k1 ||| r))))

E1 = min(e1,0.025)

e1 = t | (k2 | (s & (s1 | (k1 | r))))

E1 = min(e1,0.2)

T = min(t, 0.004)

K2 = min(k2, 0.02)

S = min(s, 0.04)

S1 = min(s1, 0.015)

K1 = min(k1, 0.02)

R = min(r, 0.02)

76/4

    19