# Routine for Multi-Stress Weibull Plotting Example
# Coded by Reuel Smith 2015-2017
# ========================================================================
# Example Problem 4.7
# Consider an accelerated test where temperature is the accelerating
# variable and 20 units were tested to failure assuming complete failures
# (no censoring). Eight units were tested at 406 K, and six units each at
# 436 K and 466 K, with times to failure tabulated below. Assuming an
# Arrhenius-Weibull life-stress relationship, find the parameters of the
# model, using both the plotting method and MLE method and compare the
# results. Find the expected life at the use level temperature of 353 K
# and compare the results.
#
# Temperature Time to Failure (hrs)
# =========== =====================================
# 406 K 248 456 528 731 813 965 972 1528
# 436 K 164 176 289 319 386 459
# 466 K 92 105 155 184 219 235
#
# Solution:
# This routine will only cover the Weibull multi-plot at the three
# different stress levels.
# ========================================================================
# PART 1: Initialize data and constants
# ========================================================================
# sorting the times to failure in ascending order
# CHANGE these values to process your own data
TTF1 <- sort(c(248,456,528,731,813,965,972,1528))
TTF2 <- sort(c(164,176,289,319,386,459))
TTF3 <- sort(c(92,105,155,184,219,235))
# set the length of the TTF vector or the number of failures
n1 <- length(TTF1)
n2 <- length(TTF2)
n3 <- length(TTF3)
# CDF or unreliability for probability plot at Stress Levels 1, 2, and 3
i1 <- c(1:n1) # uncensored data points from 1 to n1
i2 <- c(1:n2) # uncensored data points from 1 to n2
i3 <- c(1:n3) # uncensored data points from 1 to n3
# Median plotting position (i - 0.3)/(n + 0.4). Other plotting positions
# may be assumed:
# KIMBALL - (i - 0.375)/(n + 0.25)
# MIDPOINT - (i - 0.5)/n
# MEAN - i/(n + 1)
F1 <- 100*((i1-0.3)/(n1+0.4));
F2 <- 100*((i2-0.3)/(n2+0.4));
F3 <- 100*((i3-0.3)/(n3+0.4));
# Weibull scale conversion of the Failure Percent data
F1B <- log(log(1/(1-F1/100)))
F2B <- log(log(1/(1-F2/100)))
F3B <- log(log(1/(1-F3/100)))
# reliability for probability plot
R1 <- (100-F1)/100
R2 <- (100-F2)/100
R3 <- (100-F3)/100
# Form a matrix to categorize the results of the above analysis.
# Column 1 displays the TTFs, Column 2 displays the calculated Failure
# Probability, and Column 3 displays the Reliability
#TTFtable <- matrix(c(1:(3*n)), nrow = n, ncol = 3, byrow = TRUE,dimnames = list(c("1", "2", "3", "4", "5"),c("TTF", "Failure Probaility (%)", "Reliability (%)")))
#TTFtable[,1] <- TTF
#TTFtable[,2] <- F
#TTFtable[,3] <- R*100
# Scale conversion to Weibull plot scale
# Upper and lower bounds of the Percent Failure axis in percent
fc <- c(.1,99.9)
# Weibull scale conversion of the upper and lower bounds of the Percent
# Failure axis
fcB <- log(log(1/(1-fc/100)));
# Percent Failure ticks cited between the lower and upper bounds
Pticks1 <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,c(1:10),10*c(2:9),95,99,99.9)
# Weibull scale conversion of Percent Failure ticks cited between the lower
# and upper bounds
Pticks <- log(log(1/(1-Pticks1/100)))
# Percent Failure labels for the y-axis
Pticks1label <- c(0.1,0.2,0.3,"",0.5,"","","","",1,2,3,"",5,"","","","",10*c(1:9),95,99,99.9)
# ========================================================================
# PART 2: Solve for Weibull Parameters alpha and beta
# ========================================================================
# The Weibull parameters alpha and beta may be computed by assuming a
# relationship for the reliability function R(t) = exp[-(t/alpha)^beta] or
# 1/R(t) = exp[(t/alpha)^beta]. The relation translates to:
#
# ln[ln(1/R(t))] = beta*ln(t) - beta*ln(alpha)
#
# The LM function is used to fit this linear relationship and extract
# the Weibull parameters alpha and beta. See Reliability Engineering and
# Risk Analysis: A Practical Guide Section 3.3.2.2. for method.
yfit1 <- log(log(1/R1))
yfit2 <- log(log(1/R2))
yfit3 <- log(log(1/R3))
xfit1 <- log(TTF1)
xfit2 <- log(TTF2)
xfit3 <- log(TTF3)
# Calculates the slope beta and the intercept beta*ln(alpha) needed to
# extract the Weibull parameters
wblparm1 <- lm(yfit1 ~ poly(xfit1, 1, raw=TRUE))
wblparm2 <- lm(yfit2 ~ poly(xfit2, 1, raw=TRUE))
wblparm3 <- lm(yfit3 ~ poly(xfit3, 1, raw=TRUE))
beta1 <- summary(wblparm1)$coefficients[2,1]
beta2 <- summary(wblparm2)$coefficients[2,1]
beta3 <- summary(wblparm3)$coefficients[2,1]
alpha1 = exp(-summary(wblparm1)$coefficients[1,1]/beta1)
alpha2 = exp(-summary(wblparm2)$coefficients[1,1]/beta2)
alpha3 = exp(-summary(wblparm3)$coefficients[1,1]/beta3)
# CHANGE NROW if you have a different number of TTF data sets
wblresults <- matrix(c(alpha1,beta1,alpha2,beta2,alpha3,beta3), nrow = 3, ncol = 2, byrow = TRUE,dimnames = list(c("Wbl Parameters 1","Wbl Parameters 2","Wbl Parameters 3"),c("alpha", "beta")))
betamean <- mean(wblresults[,2])
# ========================================================================
# PART 3: Weibull probability plot
# ========================================================================
# Plot points at 0.1% and 99.9%
# These have a different polynomial pair than that which was used to find
# the alpha and beta
intercept1 <- summary(wblparm1)$coefficients[1,1]
intercept2 <- summary(wblparm2)$coefficients[1,1]
intercept3 <- summary(wblparm3)$coefficients[1,1]
ttfc1 <- c(exp((log(log(1./(1-0.001)))-intercept1)/beta1),exp((log(log(1./(1-0.999)))-intercept1)/beta1))
ttfc2 <- c(exp((log(log(1./(1-0.001)))-intercept2)/beta2),exp((log(log(1./(1-0.999)))-intercept2)/beta2))
ttfc3 <- c(exp((log(log(1./(1-0.001)))-intercept3)/beta3),exp((log(log(1./(1-0.999)))-intercept3)/beta3))
totttfc <- c(ttfc1,ttfc2,ttfc3)
# Computes the upper and lower bound for the TTF axis in terms of log-time
signs1 <- c(floor(log10(min(totttfc))):ceiling(log10(max(totttfc))))
logtimes1 <- 10^signs1
logspace <- c("","","","","","","","")
Pticks1X <- c(1:(9*length(logtimes1)-8))
Pticks1X[1] <- logtimes1[1]
Pticks1Xlabel <- Pticks1X
# Calculates the tick values used for the final plot
for(i2 in 1:(length(signs1)-1)){
Pticks1X[(9*i2-7):(9*(i2+1)-8)] <- logtimes1[i2]*c(2:10)
Pticks1Xlabel[(9*i2-7):(9*(i2+1)-8)] <- c("","","","","","","","",logtimes1[i2+1])
}
# The probability plot
# DEFAULT XLAB HEADING is in "Time (hours)". When using a different unit of
# time or measurement, CHANGE XLAB HEADING BELOW. Also ADD or SUBTRACT stress
# levels as needed according to the setup
# STRESS LEVEL 1 PLOT
plot(TTF1, F1B, log="x",col="blue",
xlab="Time (hours)", ylab="Percent Failure", pch=16,
xlim=c(min(logtimes1), max(logtimes1)), ylim=c(min(fcB), max(fcB)), axes=FALSE)
lines(ttfc1, fcB, col="blue")
# STRESS LEVEL 2 PLOT
points(TTF2, F2B, col="red", pch=17)
lines(ttfc2, fcB, col="red")
# STRESS LEVEL 3 PLOT
points(TTF3, F3B, col="green", pch=18)
lines(ttfc3, fcB, col="green")
axis(1, at=Pticks1X,labels=Pticks1Xlabel, las=2, cex.axis=0.7)
axis(2, at=Pticks,labels=Pticks1label, las=2, cex.axis=0.7)
#Add horizontal grid
abline(h = Pticks, lty = 2, col = "grey")
#Add vertical grid
abline(v = Pticks1X, lty = 2, col = "grey")
# The legend heading
# The first two places signify where the legend will be located. CHANGE the
# legend labels in the LEGEND="" portion. PCH changes the label
# types and LTY changes the line types.
legend(10^signs1[1], log(log(1/(1-.99))), legend=c("Stress Level 1 - 406 deg K", "Stress Level 2 - 436 deg K", "Stress Level 3 - 466 deg K"),
col=c("blue", "red", "green"), pch=c(16,17,18), cex=0.8,
text.font=4, bg='white')
wblresults
Top