9 Noncoherent system

Post date: Apr 19, 2014 12:14:17 AM

# <<The last of four parameterizations yields some negative values for probability of unavailability.  Yes, I know that's not good.>>

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

# Compute reliability of a noncoherent (nonunate) system from sparse sample data assuming independence among the components

# wipe everything from memory

rm(list=ls())

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

# infrastructural functions

many = 10000  # increase for more accuracy

constant <- function(b) if (length(b)==1) TRUE else FALSE

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

# x[i] ~ exponential(parameter), x[i] is a nonnegative integer

nextvalue.exponential <- function(x) {n <- length(x); k = sum(x); return(gammaexponential(shape=n, rate=k))}

parameter.exponential <- function(x) {n <- length(x); k = sum(x); return(gamma(shape=n, rate=k))}

gamma <- function(shape,rate=1,scale=1/rate) {rate = 1/scale; rgamma(many, shape,rate)}

gammaexponential <- function(shape,rate=1,scale=1/rate) {rate = 1/scale; rexp(many,gamma(shape,rate))}

#qgammaexponential <- function(p, shape, rate=1, scale=1/rate) rate * ((1-p)^(-1/shape) - 1)

#rgammaexponential <- function(many, shape, rate=1, scale=1/rate) return(qgammaexponential(runif(many),shape,rate,scale))

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

env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error') # serves as env, but 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]))}  # might prefer to shuffle indices rather than sample them

kn <- function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k))))  # c-box for Bernoulli parameter:   x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)

endpoints <- function(b,side, which=1:many) if (constant(b)) return(b) else if (side) return(rightside(b)[which]) else return(leftside(b)[which])

index1 <- function(mon,side) if (mon==2) return(0) else if (mon==1) return(side) else return(1-side)

index2 <- function(mon,side) if (mon==2) return(1) else if (mon==1) return(side) else return(1-side)

corners <- function(f, pr, a1=0,a2=0,a3=0,a4=0,a5=0,a6=0,a7=0,a8=0,a9=0) {

  cat('The function will be evaluated',max(2*cumprod(ifelse(pr==0,1,pr))),'x',many,'times\n')

  if (length(pr)<9) pr = c(pr,rep(0,9-length(pr)))

  f0 = rep(Inf, many)

  f1 = rep(-Inf, many)

  for (side in 0:1)

   for (i1 in index1(pr[[1]],side):index2(pr[[1]],side)) {A1 = endpoints(a1,i1)

   for (i2 in index1(pr[[2]],side):index2(pr[[2]],side)) {A2 = endpoints(a2,i2)

   for (i3 in index1(pr[[3]],side):index2(pr[[3]],side)) {A3 = endpoints(a3,i3)

   for (i4 in index1(pr[[4]],side):index2(pr[[4]],side)) {A4 = endpoints(a4,i4)

   for (i5 in index1(pr[[5]],side):index2(pr[[5]],side)) {A5 = endpoints(a5,i5)

   for (i6 in index1(pr[[6]],side):index2(pr[[6]],side)) {A6 = endpoints(a6,i6)

   for (i7 in index1(pr[[7]],side):index2(pr[[7]],side)) {A7 = endpoints(a7,i7)

   for (i8 in index1(pr[[8]],side):index2(pr[[8]],side)) {A8 = endpoints(a8,i8)

   for (i9 in index1(pr[[9]],side):index2(pr[[9]],side)) {A9 = endpoints(a9,i9)

       ff = f(A1,A2,A3,A4,A5,A6,A7,A8,A9)

       #if (side==0) f0 = pmin(f0,ff)

       #if (side==1) f1 = pmax(f1,ff)

       f0 = pmin(f0,ff)

       f1 = pmax(f1,ff)

       }}}}}}}}}

  env(f0,f1)

  }

gone = 0

decr = 0

incr = 1

both = 2

#cat('first pass, decr:',index1(decr,0), index2(decr,0),'      second pass, decr:',index1(decr,1), index2(decr,1),    '\n')

#cat('first pass, incr:',index1(incr,0), index2(incr,0),'      second pass, incr:',index1(incr,1), index2(incr,1),    '\n')

#cat('first pass, both:',index1(both,0), index2(both,0),'      second pass, both:',index1(both,1), index2(both,1),    '\n')

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 reliability of Jackson's noncoherent (nonunate) pump system

# A pump delivers liquid into a vessel monitored by a level sensor, which informs a

# level controller that actuates the pump. The reliability of this system is nonunate.  

# Below is an elaboration of the interval calculation of the system's reliability.

f <- function(a,b,c, ...) return(a * c + b - b * c)

# does this function actually compute the failure probability, or the reliability?

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

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

pump=kn(15,300)

sensor=kn(1,2000)

controller=kn(10,200)

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

# compute reliability probability

answer = corners(f, c(both,both,both), pump,sensor,controller)

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

# display results

par(mfrow=c(3,4))

plotbox(pump); title('pump')

plotbox(sensor); title('sensor')

plotbox(controller); title('controller')

plotbox(answer); title('noncoherent system')

# apparently, this result can also be effectively even more easily

sans=f(pump,sensor,controller)

plotbox(sans,col='green',new='FALSE')

# this is a surprise that I don't understand the implications of

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

# example with REAL values

# translate input from experts into parameters

q <- function(lambda,mu) return(pexp(500, lambda+mu) / (1 + mu/lambda)) # translate lambda and mu to characterizations of the component unavailability

# do the analysis

pump          =  q(1e-3, 0.1)

sensor         =  q(2e-3, 5e-2)

controller     =  q(3e-3, 1/60)

pump; sensor; controller   #    0.00990099, 0.03846154, 0.1525342

f(pump,sensor,controller)  #    0.03410508

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

# example with INTERVAL ranges from natural uncertainty from significant digits

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

# hypothetical data

makeinterval = function(a,b) env(rep(a,many),rep(b,many))

pump = makeinterval(0.003322259, 0.02912622)

sensor = makeinterval( 0.02654867, 0.05263158)

controller = makeinterval( 0.1208633, 0.1853325)

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

# compute reliability probability

answer = corners(f, c(both,both,both), pump,sensor,controller)

# this would be incorrect:

#sans=f(pump,sensor,controller)

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

# display results

#par(mfrow=c(2,2))

plotbox(pump); title('pump')

plotbox(sensor); title('sensor')

plotbox(controller); title('controller')

plotbox(answer); title('noncoherent system')

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

# interval outputs agree with results computed by Risk Calc

#pump          = q1 = q([1e-3 per hour], [0.1 per hour])

#sensor         = q2 = q([2e-3 per hour], [5e-2 per hour])

#controller     = q3 = q([3e-3 per hour], 1/[60 hours])       +   0

#

#[1e-3 per hour]; [0.1 per hour]

#    [ 0.0005, 0.0015] per hour

#    [ 0.05, 0.15] per hour

#[2e-3 per hour]; [5e-2 per hour]

#    [ 0.0015, 0.0025] per hour

#    [ 0.045, 0.055] per hour

#[3e-3 per hour]; [60 hours]

#    [ 0.0025, 0.0035] per hour

#    [ 55, 65] hours

#

# pump; sensor; controller

#    [ 0.003322259, 0.02912622]

#    [ 0.02654867, 0.05263158]

#    [ 0.1208633, 0.1853325]

# answer

#    [ 0.02224406, 0.04979065]

#

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

# hypothetical data in the form of failure times informing exponential C-BOXES then translated into unavailability

# parameter.exponential(data) gives c-box for lambda of exponential distribution based on empirical failure times

# nextvalue.exponential(data) gives p-box of future exponential failure times

pdata = rexp(7, 1e-3)

sdata = rexp(17, 2e-3)

cdata = rexp(71, 3e-3)

pump = nextvalue.exponential(pdata)

sensor = nextvalue.exponential(sdata)

controller = nextvalue.exponential(cdata)

# how should we translate lambda and mu to characterizations of the component unavailability?

q <- function(lambda,mu) return(pexp(500, lambda+mu) / (1 + mu/lambda))

#pump          = q1 = q([1e-3 per hour], [0.1 per hour])

#sensor         = q2 = q([2e-3 per hour], [5e-2 per hour])

#controller     = q3 = q([3e-3 per hour], 1/[60 hours])       +   0

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

# compute reliability probability

answer = corners(f, c(both,both,both), pump,sensor,controller)

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

# display results

#par(mfrow=c(3,4))

plotbox(pump); title('pump')

plotbox(sensor); title('sensor')

plotbox(controller); title('controller')

plotbox(answer); title('noncoherent system')

lines(c(0,0),c(0,1))