# R & WinBUGS Routine to Perform a Step-Stress Testing Example
# Coded by Reuel Smith 2015-2017
#
# NOTE:
# Install and/or Load the R2WINBUGS Package prior to using this routine
# (Or R2OPENBUGS depending on what program you use. WinBUGS or OpenBUGS
# must be installed on the computer prior to using this routine.
# ========================================================================
# Example Problem 4.19
# Consider a step-stress test of cable insulation. The test was run to
# estimate life at a design stress of 400 volts/mm. The data in Table 4.24
# illustrate the data analysis methods. Each specimen was first held for
# 10 minutes each at 5 kV, 10 kV, 15 kV, and 20 kV before it went into step
# 1 at 26 kV. In steps 1 through 10, a specimen has the same hold time at
# each voltage (15 minutes, 1 hour, 4 hours, or 16 hours). Table 4.24
# shows the step number and the total time on test when a specimen broke
# down, and its insulation thickness (used to calculate its stress as the
# voltage divided by the thickness). Assume that this tested cable 1 is
# comparable to another cable 2 tested before, so we can use the MLE result
# of the test of cable 2 as the prior information. Assume that cable 2 was
# tested and the data shown in Table 4.25 were obtained. Update the models
# parameters using the Bayesian estimation.
#
# Solution:
#
# ========================================================================
# Set the working directory
setwd("C:/Users/ReuelS/Documents/ENRE447")
# ========================================================================
# Initialize data and constants
# ========================================================================
# Definitions of different variables used in this code
# x(1): A
# x(2): p
# x(3): beta
# Voltage pattern in volts by step number
V <- 1000*c(26,28.5,31,33.4,36,38.5,41,43.5,46,48.5)
# Total number of steps
STEPS <- 10
# Cable data:
# First Column - The hold time for each step. If the final step is 5, then
# the total time elapsed will be Hold Time x 5
# Second Column - Final Step which is the step that testing ends either
# because of specimen failure or for censoring
# Third Column - Total Time-to-Failure. Times-to-failure that are given as
# negative are treated as censored data or a runtime without failure.
# Fourth Column - Thickness of the cable insulation (mm).
# Cable data for first one (will be for updating of priors)
cable1 <- matrix(c(15,5,102,27,0,0,0,0,
15,5,113,27,0,0,0,0,
15,5,113,27,0,0,0,0,
60,6,-370,29.5,0,0,0,0,
60,6,345,29.5,0,0,0,0,
60,6,-345,28,0,0,0,0,
240,6,1249,29,0,0,0,0,
240,6,1333,29,0,0,0,0,
240,6,-1333,29,0,0,0,0,
240,5,1096.9,29,0,0,0,0,
240,6,1250.8,30,0,0,0,0,
240,5,1097.9,29,0,0,0,0,
960,3,2460.9,30,0,0,0,0,
960,3,-2460.9,30,0,0,0,0,
960,3,2700.4,30,0,0,0,0,
960,4,2923.9,30,0,0,0,0,
960,2,1160,30,0,0,0,0,
960,3,1962.9,30,0,0,0,0,
960,1,-363.9,30,0,0,0,0,
960,1,-898.4,30,0,0,0,0,
960,5,4142.1,30,0,0,0,0), nrow = 21, ncol = 8, byrow = TRUE)
# Cable data for second one (will be used to obtain the priors
cable2 <- matrix(c(15,5,114,27,0,0,0,0,
15,6,119,27,0,0,0,0,
15,6,116,27,0,0,0,0,
15,6,117,27,0,0,0,0,
60,7,423,29.5,0,0,0,0,
60,7,-423,28,0,0,0,0,
240,6,1401.5,29,0,0,0,0,
240,6,1261.7,29,0,0,0,0,
240,6,1419.4,29,0,0,0,0,
240,6,1364.2,29,0,0,0,0,
960,2,1398.1,30,0,0,0,0,
960,3,2691.7,30,0,0,0,0,
960,3,2000,30,0,0,0,0,
960,3,-2000,30,0,0,0,0,
960,3,2149.6,30,0,0,0,0,
960,2,1429.7,30,0,0,0,0,
960,2,1502.7,30,0,0,0,0,
960,2,1440.9,30,0,0,0,0,
960,4,3349.1,30,0,0,0,0,
960,4,-3349.1,30,0,0,0,0), nrow = 20, ncol = 8, byrow = TRUE)
# The initial time in minutes. Before step 1, each specimen was held 10
# minutes each at 5, 10, 15, and 20 kV, so 40 minutes passed before the
# start of the step-stress testing.
t0 = 40
# ========================================================================
# Separation of Failure and Survivor Data for Cable 1
# ========================================================================
# This loop takes the existing cable data and computes the stress levels of
# each data (in V/mm) and the true start time based on the initial testing
# time of 40 minutes. The true start time is computed from based on the
# final step n t_(n-1)
cable1[,6]<-t0 + (STEPS-1)*cable1[,1]
cable1[,8]<-V[10]/cable1[,4]
for(i in 1:length(cable1[,1])) {
# Stress levels of Cable 1 in V/mm (i)
cable1[i,7]<-V[cable1[i,2]]/cable1[i,4]
# Previous time on Cable 1 in minutes (i-1)
if(cable1[i,2]>1) {
cable1[i,5]=t0+(cable1[i,2]-1)*cable1[i,1]
}
else {
cable1[i,5]=t0
}
}
# This loop separates the failure data from the survivor data for
# preparation of the Bayesian analysis step
i_F<-1
i_S<-1
fail1 <- 0*c(1:length(cable1[,1]))
surv1 <- 0*c(1:length(cable1[,1]))
for(i in 1:length(cable1[,1])){
if(cable1[i,3]>0){
fail1[i] <- 1
i_F <- i_F + 1
}
else{
surv1[i] <- 1
i_S <- i_S + 1
}
}
cable1F <- matrix(c(1:(i_F-1)*8), nrow = (i_F-1), ncol = 8, byrow = TRUE)
cable1S <- matrix(c(1:(i_S-1)*8), nrow = (i_S-1), ncol = 8, byrow = TRUE)
i_F<-1
i_S<-1
for(i in 1:length(cable1[,1])){
if(cable1[i,3]>0){
cable1F[i_F,] <- cable1[i,]
i_F <- i_F + 1
}
else{
cable1S[i_S,] <- abs(cable1[i,])
i_S <- i_S + 1
}
}
# ========================================================================
# Separation of Failure and Survivor Data for Cable 2
# ========================================================================
# This loop takes the existing cable data and computes the stress levels of
# each data (in V/mm) and the true start time based on the initial testing
# time of 40 minutes. The true start time is computed from based on the
# final step n t_(n-1)
cable2[,6]<-t0 + (STEPS-1)*cable2[,1]
cable2[,8]<-V[10]/cable2[,4]
for(i in 1:length(cable2[,1])) {
# Stress levels of Cable 2 in V/mm (i)
cable2[i,7]<-V[cable2[i,2]]/cable2[i,4]
# Previous time on Cable 2 in minutes (i-1)
if(cable2[i,2]>1) {
cable2[i,5]=t0+(cable2[i,2]-1)*cable2[i,1]
}
else {
cable2[i,5]=t0
}
}
# This loop separates the failure data from the survivor data for
# preparation of the Bayesian analysis step
i_F<-1
i_S<-1
fail2 <- 0*c(1:length(cable2[,1]))
surv2 <- 0*c(1:length(cable2[,1]))
for(i in 1:length(cable2[,1])){
if(cable2[i,3]>0){
fail2[i] <- 1
i_F <- i_F + 1
}
else{
surv2[i] <- 1
i_S <- i_S + 1
}
}
cable2F <- matrix(c(1:(i_F-1)*8), nrow = (i_F-1), ncol = 8, byrow = TRUE)
cable2S <- matrix(c(1:(i_S-1)*8), nrow = (i_S-1), ncol = 8, byrow = TRUE)
i_F<-1
i_S<-1
for(i in 1:length(cable2[,1])){
if(cable2[i,3]>0){
cable2F[i_F,] <- cable2[i,]
i_F <- i_F + 1
}
else{
cable2S[i_S,] <- abs(cable2[i,])
i_S <- i_S + 1
}
}
# ========================================================================
# Initialize Formulae for evaluation
# ========================================================================
# SETUP FOR TAU (tau1 through tau9 for stress levels 2 through 10. The
# tau values calculated follow the equation:
#
# tau_(i-1)=(t_(i-1)-t_(i-2)+tau_(i-2) ) (V_(i-1)/V_i )^p
#
# The (V(i-1)/V(i))^p term is defined as the inverse of the acceleration
# factor (AF) between stress V(i-1) and the next stress V(i) under the
# Inverse Power Relation (IPR) for life (or 1/AF = 1/[(V(i)/V(i-1))^p]).
# This value will have to be changed in the tau statements when and if
# assuming a different life model. Parameters being sought must all be
# given an x(*) value:
# Example: IPR - x(1)=A, x(2)=p, x(3)=beta
# ========================================================================
# Acceleration factors: Change (V(i-1)/V(i))^p statement in tau for
# different models:
#
# Inverse Power Law: (V(i-1)/V(i))^p
# Eqn(4.69), V-stress, p-parameter
#
# Arrhenius: exp(a*(1/T(i) - 1/T(i-1)))
# Eqn(4.40), T-temperature, a-slope
#
# Eyring: (V(i-1)/V(i))*exp(a*(1/T(i) - 1/T(i-1)))
# Eqn(4.55), V-stress, T-temperature, a-slope
#
# Dual Stress: exp(a*(1/T(i) - 1/T(i-1)) + b*(1/V(i) - 1/V(i-1)))
# Eqn(4.85), T-temperature, V-non-thermal stress, a-slope 1, b-slope 2
#
# Power Exponential: ((V(i-1)/V(i))^n)*exp(a*(1/T(i) - 1/T(i-1)))
# Eqn(4.97), V-stress, T-temperature, a-slope 1, n-slope 2
# ========================================================================
# The tau values for all 10 steps
# The input for each tau function is the parameters (for parameter p) and
# the hold time for each step (n1). Taus are calculated up to the final
# step STEPS = 10.
# ========================================================================
# tau for step 2
tau1 <- function(n,S) {
(t0 + 1*S)*(V[1]/V[2])^n
}
# tau for step 3
tau2 <- function(n,S) {
((t0 + 2*S)-(t0 + 1*S)+tau1(n,S))*(V[2]/V[3])^n
}
# tau for step 4
tau3 <- function(n,S) {
((t0 + 3*S)-(t0 + 2*S)+tau2(n,S))*(V[3]/V[4])^n
}
# tau for step 5
tau4 <- function(n,S) {
((t0 + 4*S)-(t0 + 3*S)+tau3(n,S))*(V[4]/V[5])^n
}
# tau for step 6
tau5 <- function(n,S) {
((t0 + 5*S)-(t0 + 4*S)+tau4(n,S))*(V[5]/V[6])^n
}
# tau for step 7
tau6 <- function(n,S) {
((t0 + 6*S)-(t0 + 5*S)+tau5(n,S))*(V[6]/V[7])^n
}
# tau for step 8
tau7 <- function(n,S) {
((t0 + 7*S)-(t0 + 6*S)+tau6(n,S))*(V[7]/V[8])^n
}
# tau for step 9
tau8 <- function(n,S) {
((t0 + 8*S)-(t0 + 7*S)+tau7(n,S))*(V[8]/V[9])^n
}
# tau for step 10
tau9 <- function(n,S) {
((t0 + 9*S)-(t0 + 8*S)+tau8(n,S))*(V[9]/V[10])^n
}
# The tau values and adjusted times for for Cable 1. The tau functions
# that go into tauC1 must be predetermined according to the step number
# (for failure times only). This is always one number less than the given
# step. For example: at step 2 use tau 1; for step 5 use tau4.
# Adjusted time for failures
adjtimeC1F <- function(n) {
(cable1F[,3] - cable1F[,5])*((cable1F[,8]/cable1F[,7])^n)
}
adjtimeC2F <- function(n) {
(cable2F[,3] - cable2F[,5])*((cable2F[,8]/cable2F[,7])^n)
}
# Adjusted time for survivals
adjtimeC1S <- function(n) {
(cable1S[,3] - cable1S[,5])*((cable1S[,8]/cable1S[,7])^n)
}
adjtimeC2S <- function(n) {
(cable2S[,3] - cable2S[,5])*((cable2S[,8]/cable2S[,7])^n)
}
# Adjusted time for Cable 1
adjtimeC1 <- function(n) {
(abs(cable1[,3]) - cable1[,5])*((cable1[,8]/cable1[,7])^n)
}
# Adjusted time for Cable 2
adjtimeC2 <- function(n) {
(abs(cable2[,3]) - cable2[,5])*((cable2[,8]/cable2[,7])^n)
}
# ====================================================================
# SETUP FOR SAMPLING - CABLE 2
# ====================================================================
niter <- 600000
# Parameter 1 Priors: A_logsig and A_tau. Parameter 1 A is assumed to have
# a lognormal prior distribution LOGN(A_logsig,A_tau). A_tau is
# 1/A_logsig^2 as per the definition given in WinBUGS
A_mu <- 0.00000000001
A_sig <- 100
A_logmu <- log(A_mu/sqrt(1 + ((A_sig/A_mu)^2)))
A_logsig <-sqrt(log(1 + ((A_sig/A_mu)^2)))
A_tau <- 1/(A_logsig^2)
# Parameter 2 Priors: p_L and p_H. Parameter 2 p is assumed to have a
# uniform prior distribution UNIF(p_L,p_H).
p_L <- -1000
p_H <- 1000
# Parameter 3 Priors: beta_L and beta_H. Parameter 2 p is assumed to have a
# uniform prior distribution UNIF(beta_L,beta_H).
beta_L <- 0
beta_H <- 1000
t <- abs(cable2[,3])
tS <- cable2[,5]
nt <- length(t)
Si <- cable2[,7]
S <- cable2[,8]
fail <- fail2
surv <- surv2
data <- list("A_logmu","A_tau","p_L","p_H","beta_L","beta_H","t","tS","nt","S","Si","fail","surv")
parameters <- c("A","p","beta")
inits <- list(list(A = 0.00000000001, p = 1, beta = 1))
# Under MODEL.FILE, change the folder name to where your WinBUGS model
# is located. Under BUG.DIRECTORY, change the folder name to where the
# WinBUGS executable is located.
ADT.sim <- bugs(data,inits,parameters,model.file="C:/Users/ReuelS/Documents/ENRE447/StepStressExample_4_19.odc",n.chains = 1, n.iter = niter, n.burnin=25000, n.thin=16, debug=TRUE, bugs.directory="C:/Users/ReuelS/Desktop/POD Research/WinBUGS/winbugs14/WinBUGS14", DIC=FALSE)
attach.bugs(ADT.sim)
#detach.bugs()
flush.console()
A2<-A
p2<-p
beta2<-beta
results1 <- matrix(c(min(A2),median(A2),max(A2),mean(A2),sd(A2),min(p2),median(p2),max(p2),mean(p2),sd(p2),min(beta2),median(beta2),max(beta2),mean(beta2),sd(beta2)), nrow = 3, ncol = 5, byrow = TRUE,dimnames = list(c("A", "p", "beta"),c("min", "median", "max", "mean", "SD")))
# ====================================================================
# SETUP FOR SAMPLING - CABLE 1
# ====================================================================
niter <- 600000
# Parameter 1 Priors: A_logsig and A_tau. Parameter 1 A is assumed to have
# a lognormal prior distribution LOGN(A_logsig,A_tau). A_tau is
# 1/A_logsig^2 as per the definition given in WinBUGS
A_mu <- mean(A2)
A_sig <- sd(A2)
A_logmu <- log(A_mu/sqrt(1 + ((A_sig/A_mu)^2)))
A_logsig <-sqrt(log(1 + ((A_sig/A_mu)^2)))
A_tau <- 1/(A_logsig^2)
# Parameter 2 Priors: p_L and p_H. Parameter 2 p is assumed to have a
# uniform prior distribution UNIF(p_L,p_H).
p_L <- min(p2)
p_H <- max(p2)
# Parameter 3 Priors: beta_L and beta_H. Parameter 2 p is assumed to have a
# uniform prior distribution UNIF(beta_L,beta_H).
beta_L <- min(beta)
beta_H <- max(beta)
t <- abs(cable1[,3])
tS <- cable1[,5]
nt <- length(t)
Si <- cable1[,7]
S <- cable1[,8]
fail <- fail1
surv <- surv1
data <- list("A_logmu","A_tau","p_L","p_H","beta_L","beta_H","t","tS","nt","S","Si","fail","surv")
parameters <- c("A","p","beta")
inits <- list(list(A <- rlnorm(1,A_logmu,A_logsig), p <- runif(1,p_L,p_H), beta <- runif(1,beta_L,beta_H)))
# Under MODEL.FILE, change the folder name to where your WinBUGS model
# is located. Under BUG.DIRECTORY, change the folder name to where the
# WinBUGS executable is located.
ADT.sim <- bugs(data,inits,parameters,model.file="C:/Users/ReuelS/Documents/ENRE447/StepStressExample_4_19.odc",n.chains = 1, n.iter = niter, n.burnin=25000, n.thin=16, debug=TRUE, bugs.directory="C:/Users/ReuelS/Desktop/POD Research/WinBUGS/winbugs14/WinBUGS14", DIC=FALSE)
attach.bugs(ADT.sim)
#detach.bugs()
flush.console()
A1<-A
p1<-p
beta1<-beta
results2 <- matrix(c(min(A1),median(A1),max(A1),mean(A1),sd(A1),min(p1),median(p1),max(p1),mean(p1),sd(p1),min(beta1),median(beta1),max(beta1),mean(beta1),sd(beta1)), nrow = 3, ncol = 5, byrow = TRUE,dimnames = list(c("A", "p", "beta"),c("min", "median", "max", "mean", "SD")))
results1
results2
Top