Bullmer Fabric Cutter Reliability

library(WeibullR)
library(readxl)
rm(list = ls()) # remove all environment
Blmr <- read_excel("Bullmer.xlsx", col_types = c("numeric", "skip", "text", "text", "text", "text", 
                                                    "numeric", "numeric", "numeric", "skip", 
                                                    "skip", "skip"))
#View(Blmr)

f <- function(time, title, pch) {
  today = 3050
  
  d <- wblr(time, pp = "median", pch = pch)
  d$options$main <- c(title , d$options$main)
  d2 <- wblr.conf(wblr.fit(d, dist = "weibull2p"))
  d3 <- wblr.fit(d, dist = "weibull3p")
  e2 <- wblr.conf(wblr.fit(d, dist = "lognormal"))
  e3 <- wblr.fit(d, dist = "lognormal3p")
  
  time<-c(sort(time),today)
  w <- c(d2$fit[[1]]$beta,d2$fit[[1]]$eta)
  haz2 <- w["Beta"]/w["Eta"]*(time/w["Eta"])^(w["Beta"]-1)
  w <- c(d3$fit[[1]]$beta,d3$fit[[1]]$eta, d3$fit[[1]]$t0)
  haz3 <- w["Beta"]/w["Eta"]*((time-w["t0"])/w["Eta"])^(w["Beta"]-1)
  
  w4 <- e2$fit[[1]]$fit_vec
  w5 <- e3$fit[[1]]$fit_vec
  haz4 <- dlnorm(time, w4["Mulog"], w4["Sigmalog"],log = F)/(plnorm(time,w4["Mulog"], w4["Sigmalog"],lower.tail = F, log.p = F))
  haz5 <- dlnorm(time, w5["Mulog"], w5["Sigmalog"],log = F)/(plnorm(time,w5["Mulog"], w5["Sigmalog"],lower.tail = F, log.p = F))
  
  LEGEND<-c(paste("Weibull 2p:",format(d2$fit[[1]]$fit_vec[3],digits = 3)), 
            paste("Weibull 3p:",format(d3$fit[[1]]$fit_vec[4],digits = 3)),
            paste("Lognormal 2p:",format(w4[3], digits = 3)),
            paste("Lognormal 3p:",format(w5[4], digits = 3)))
  plot.wblr(d2, in.legend.blives = F)
  plot.wblr(d3)
  plot.wblr(e2, in.legend.blives = F)
  plot.wblr(e3)
  
  if (w["t0"]>0) {
    YLIM<-c(min(haz3,haz2,haz4,haz5),max(haz3,haz2,haz4,haz5))
  } else { YLIM<-c(min(haz2,haz4),max(haz2,haz4)) }
  YLIM<-YLIM*100
  plot(time,haz2*100, type = "b", main = c(title,"Hazard rate"), xlab = "Time", ylab="Failure rate [%]",col="blue", ylim = YLIM, pch=0)
  lines(time,haz4*100, type = "b", col="green", lty=3,pch=2)
  if (w["t0"]>0) {
    lines(time,haz3*100, type = "b", col="red", lty=2,pch=1 ) 
    lines(time,haz5*100, type = "b" ,col="black" , lty=4,pch=3)
  }
  grid()
  legend(x="topleft",legend = LEGEND, cex = 0.7, col=c("blue","red","green","black"),pch=c(0,1,2,3))
}

f2<-function(frame, title) {
  mldt<-mean(frame$LDT , na.rm = T)
  mttr<-mean(frame$TTR , na.rm = T)
  mtbf<-f3(frame$Event)
  print(paste(title," Availability:"))
  print(paste("MTBF:",mtbf))
  print(paste("MTTR:",mttr))
  print(paste("MLDT:",mldt))
  print(paste("Technical:",format(mtbf/(mtbf+mttr),digits = 3),"[%]"))
  print(paste("Operational:",format(mtbf/(mtbf+mldt+mttr),digits = 3),"[%]"))
}

f3<-function(frame) {
  frame<-sort(frame)
  l<-length(frame)
  x<-0.0
  for (i in 2:l) {
    x<-x+frame[i]-frame[i-1]
  }
  return(x/(l-1))
}
q <- wblr(Blmr$Event, pch = 0)
q1 <- wblr.conf(wblr.fit(q, dist = "weibull2p"))
plot.wblr(q1, in.legend.blives = F)

q2 <- wblr.conf(wblr.fit(q, dist = "lognormal"))
plot.wblr(q2, in.legend.blives = F)

q3 <- wblr.fit(q, dist = "weibull3p")
plot.wblr(q3)

q4 <- wblr.fit(q, dist = "lognormal3p")
plot.wblr(q4)

f2(Blmr,"Bulmer fabric cutter")
## [1] "Bulmer fabric cutter  Availability:"
## [1] "MTBF: 35.0135135135135"
## [1] "MTTR: 2.12500000000029"
## [1] "MLDT: 1.96"
## [1] "Technical: 0.943 [%]"
## [1] "Operational: 0.896 [%]"
f(time = Blmr$Event[28:38], title = "Control events", pch = 2)

f2(Blmr[28:38,],"Control")
## [1] "Control  Availability:"
## [1] "MTBF: 30"
## [1] "MTTR: 0.0916666666671517"
## [1] "MLDT: 0"
## [1] "Technical: 0.997 [%]"
## [1] "Operational: 0.997 [%]"
f(time = Blmr$Event[39:42], title = "Cutter failure events", pch = 3)

f2(Blmr[39:43,],"Cutter")
## [1] "Cutter  Availability:"
## [1] "MTBF: 600.5"
## [1] "MTTR: 0.71875"
## [1] "MLDT: 0"
## [1] "Technical: 0.999 [%]"
## [1] "Operational: 0.999 [%]"
f(time = Blmr$Event[44:50], title = "General Electric events", pch = 4)

f2(Blmr[44:50,],"Electric")
## [1] "Electric  Availability:"
## [1] "MTBF: 379.166666666667"
## [1] "MTTR: 0.0625"
## [1] "MLDT: 0"
## [1] "Technical: 1 [%]"
## [1] "Operational: 1 [%]"
f(time = Blmr$Event[51:55], pch = 5, title = "Head System events")

f2(Blmr[51:55,],"Head")
## [1] "Head  Availability:"
## [1] "MTBF: 565"
## [1] "MTTR: 0.145833333335759"
## [1] "MLDT: 0"
## [1] "Technical: 1 [%]"
## [1] "Operational: 1 [%]"
f(time = Blmr$Event[57:63], pch = 6, title = "Software related events")

f2(Blmr[57:63,],"Software")
## [1] "Software  Availability:"
## [1] "MTBF: 291.833333333333"
## [1] "MTTR: 0.527777777778586"
## [1] "MLDT: 0.333333333333333"
## [1] "Technical: 0.998 [%]"
## [1] "Operational: 0.997 [%]"
f(time = Blmr$Event[64:67], pch = 8, title = "Stopper related events")

f2(Blmr[64:67,],"Stopper")
## [1] "Stopper  Availability:"
## [1] "MTBF: 258.333333333333"
## [1] "MTTR: 10.1805555555523"
## [1] "MLDT: 10"
## [1] "Technical: 0.962 [%]"
## [1] "Operational: 0.928 [%]"
f(time = Blmr$Event[68:75], pch = 9, title = "Vacuum System events")

f2(Blmr[68:75,],"Vacuum")
## [1] "Vacuum  Availability:"
## [1] "MTBF: 302.285714285714"
## [1] "MTTR: NaN"
## [1] "MLDT: NaN"
## [1] "Technical: NaN [%]"
## [1] "Operational: NaN [%]"