---
title: "Reliability Testing"
output: 
  learnr::tutorial:
    progressive: true
    allow_skip: true
runtime: shiny_prerendered
---

```{r setup, include=FALSE}
library(learnr)
```

## Introduction

Welcome to the tutorial on **Reliability Testing**! In this tutorial, you will learn about different types of reliability tests and how to apply them. These models are commonly applied in industries like aerospace, automotive, and electronics, where improving product reliability over time is critical.

## What is Reliability Testing?

Reliability testing is a process to ensure that a product or system performs its intended functions under specified conditions over a defined period. For example, in the automotive industry, we may want to design a car that runs without breaking down for 5 to 10 years. 

## Why is Reliability Testing important?

Reliability testing helps organizations predict the future reliability of products, reduce failure rates, and lower costs associated with repairs and downtime. 

## Learning Objectives

By the end of this tutorial, you will be able to:

1. Define key reliability growth concepts, including Crow-AMSAA and Duane models.
2. Fit a reliability growth model to real-world data using R.
3. Interpret reliability growth plots and identify trends.
4. Apply the Crow-AMSAA model to assess reliability growth.
5. Explain fundamental concepts of accelerated life testing, including the use of Arrhenius and Power Law Models.
6. Conduct an accelerated life test with real-world datasets, utilizing R for analysis.
7. Analyze plots that illustrate the relationships in accelerated life testing, identifying key patterns and data trends.
8. Utilize Arrhenius and Power Law models to evaluate the impact of stress factors on product reliability.

## Reliability Growth Analysis

**Reliability Growth Analysis (RGA)** is a method to monitor and improve, or grow, the reliability of a product or system over time. Often the product is improved due to changes in the design throughout the testing period.

## The Duane Model

**The Duane Model** was developed by J. T. Duane in 1964 and provides a simple and graphical way to observe and analyze whether failure rates are improving as changes are made to a product or system. 

Specifically, the Duane Model is a log-log plot of the cumulative Mean Time Between Failures (MTBF) vs cumulative time.

The slope of the line on the plot indicates the rate of reliability growth:

* A positive slope means that the system is improving (reliability is growing, the failure rate is decreasing).
* A zero slope means there is no change in reliability (the system is stable).
* A negative slope indicates that reliability is worsening (the failure rate is increasing).

### The Duane Model: Reliability Growth vs Deteriation

Take the following example Duane Plot. Adjust the Beta value input to see the failure rate both decreasing and increasing.

```{r echo=FALSE}
numericInput("beta", label = h3("Beta Value:"), value = 1, min = 0, max = 5, step = 0.1, width = '50%')
plotOutput("duanePlot")
```

```{r, context="server"}
output$duanePlot <- renderPlot({
  x <- seq(1, 100, by = .1)
  y <- 0.1 * x ^ (input$beta - 1)
  log_x <- log(x)
  log_y <- log(y)
  mainTitle <- if (input$beta == 1) {
    "Duane Plot with Stable Failure Rate"
  } else if (input$beta < 1) {
    "Duane Plot with Increasing Failure Rate"
  } else {
    "Duane Plot with Decreasing Failure Rate"
  }
  plot(log_x, log_y, type = "l", col = "blue", xlab = "Time", ylab = "Cumulative MTBF", main = mainTitle)
})
```

### The Duane Model Equation

The cumulative MTBF at any time \( t \) can be expressed as:

\[
CMBTF(t) = K \cdot t^{\beta - 1}
\]

Where:

- \( CMBTF(t) \) is the cumulative MTBF at time \( t \),
- \( K \) is a constant related to the initial failure rate,
- \( t \) is the cumulative operating time,
- \( \beta \) is the reliability growth parameter (slope of the line on the log-log plot).

  - If \( \beta > 1 \), the system reliability is improving (the failure rate decreases over time).
  - If \( \beta = 1 \), the system reliability is constant (no improvement or degradation).
  - If \( \beta < 1 \), the system reliability is deteriorating (the failure rate increases over time).

### Exercises

```{r duane1, echo=FALSE}
question("What does a Duane Plot with a slope of greater than 1 indicate?",
  answer("A) The system's reliability is improving", correct = TRUE),
  answer("B) The system's reliability is deteriorating"),
  answer("C) The system has no change in reliability"),
  random_answer_order = TRUE
)
```

```{r duane2, echo=FALSE}
question("Which beta value may indicate an infant mortality failure mode?",
  answer("1"),
  answer("0.5"),
  answer("2"),
  answer("Unable to be determined in this case", correct = TRUE),
  random_answer_order = TRUE
)
```

## ReliaGrowR

### Getting Started with ReliaGrowR

For this tutorial, we will use **ReliaGrowR**: a R package for Reliability Growth Analysis.

First, check if ReliaGrowR is installed in R and install if not.

```{r echo=TRUE, message=FALSE, warning=FALSE}
if(!(require(ReliaGrowR))){install.packages('ReliaGrowR')}
```

Next, set up some cumulative time and failure data. Then use the `duane_plot` function to fit a Duane Model to the data. 

```{r echo=TRUE}
times <- c(100, 200, 300, 400, 500)
failures <- c(1, 2, 1, 3, 2)
fit <- duane_plot(times, failures, line_col = "red")
```

### Exercises

```{r duane3, echo=FALSE}
question("In this example, what is the approximate cumulative MTBF at time 1000?",
  answer("200"),
  answer("100"),
  answer("120"),
  answer("150", correct = TRUE),
  random_answer_order = TRUE
)
```

```{r duane4, echo=FALSE}
question("How would you describe the reliability of this system?",
  answer("The system reliability is improving", correct = TRUE),
  answer("The system reliability is worsening"),
  answer("The system reliability is Stable"),
  random_answer_order = TRUE
)
```

## The Crow-AMSAA Model

**The Crow-AMSAA** is a **Non-Homogeneous Poisson Process (NHPP)** model developed by Larry H. Crow and is based on the Army Materiel Systems Analysis Activity (AMSAA) testing framework.

Specifically, the Crow-AMSAA Model is a plot of the cumulative number of failures vs cumulative time. 

Similar to the Duane Model, the shape of the model indicates the rate of reliability growth:

* A model fit with downward curvature means that the system is improving (reliability is growing, the failure rate is decreasing).
* A linear model means there is no change in reliability (the system is stable).
* A model fit with upward curvature indicates that reliability is worsening (the failure rate is increasing).

### The Crow-AMSAA Model: Reliability Growth vs Deteriation

Take the following example Crow-AMSAA model. Adjust the Beta value input to see the failure rate both decreasing and increasing.

```{r echo=FALSE}
numericInput("beta", label = h3("Beta Value:"), value = 1, min = 0, max = 5, step = 0.1, width = '50%')
plotOutput("crowPlot")
```

```{r, context="server"}
output$crowPlot <- renderPlot({
  x <- seq(1, 100, by = .1)
  y <- 0.5 * x ^ input$beta
  mainTitle <- if (input$beta == 1) {
    "Crow-AMSAA Model with Stable Failure Rate"
  } else if (input$beta < 1) {
    "Crow-AMSAA Model with Decreasing Failure Rate"
  } else {
    "Crow-AMSAA Model with Increasing Failure Rate"
  }
  plot(x, y, type = "l", col = "darkgreen", xlab = "Time", ylab = "Cumulative Failures", main = mainTitle)
})
```

### The Crow-AMSAA Model Equation:

The cumulative number of failures up to time \( t \) is given by:

\[N(t) = \lambda_0 \cdot t^{\beta}\]

Where:

* \( N(t) \) is the cumulative number of failures,
* \( \lambda_0 \) is the scale parameter.

* The Shape Parameter (\( \beta \)) indicates whether the system is improving or deteriorating:

  * If \( \beta \) > 1, failures are increasing over time, indicating that reliability is worsening.
  * If \( \beta \) < 1, failures are decreasing, indicating that reliability is improving over time.
  * \( \beta \) = 1 implies a constant failure rate (no growth or degradation).

### Crow-AMSAA with ReliaGrowR

To run a Crow-AMSAA model with ReliaGrowR, first set up some cumulative time and failure data:
```{r}
times <- c(100, 200, 300, 400, 500)
failures <- c(1, 2, 1, 3, 2)
```

Then use the `rga` function to fit the model and the `plot_rga` function to plot the results:
```{r}
result <- rga(times, failures)
plot_rga(result, line_col = "orange")
```

### Exercises

```{r rga1, echo=FALSE}
question("In this example, what is the cumulative number of failures at time 400?",
  answer("4"),
  answer("8"),
  answer("6"),
  answer("7", correct = TRUE),
  random_answer_order = TRUE
)
```

```{r rga2, echo=FALSE}
question("How would you describe the reliability of this system?",
  answer("The reliability of the system is improving", correct = TRUE),
  answer("The reliability of the system is worsening"),
  answer("The reliability of the system is stable"),
  random_answer_order = TRUE
)
```

### The Piecewise NHPP Model

The **Piecewise NHPP** model is an extension of the standard NHPP model that includes different segments or phases of time that follow separate failure distributions. This model is particularly useful when a system experiences changes in failure behavior over different development phases, such as the initial, interim and final phases of a development process.

### The Failure Distribution

The failure intensity function for each time segment \(i\) is modeled as:

\[\lambda_i(t) = \frac{\beta_i}{\eta_i} \left( \frac{t - t_{i-1}}{\eta_i} \right)^{\beta_i - 1}\]

Where:

* \( \lambda_i(t) \) is the failure intensity in time interval \( i \),
* \( \beta_i \) is the shape parameter for interval \( i \),
* \( \eta_i \) is the scale parameter (characteristic life) for interval \( i \),
* \( t_{i-1} \) is the start time of the interval.

### Example

To run a Piecewise NHPP with ReliaGrowR, first, set up some cumulative time and failure data and specify a breakpoint:
```{r}
times <- c(25, 55, 97, 146, 201, 268, 341, 423, 513, 609, 710, 820, 940, 1072, 1217)
failures <- c(1, 1, 2, 4, 4, 1, 1, 2, 1, 4, 1, 1, 3, 3, 4)
breaks <- 500
```

Then run the rga and plot the results:
```{r}
result <- rga(times, failures, model_type = "Piecewise NHPP", breaks = breaks)
plot_rga(result)
```

### The Piecewise NHPP with Change Point Detection

The **Piecewise NHPP with Change Point Detection** is an advanced model to identify changes in failure behavior and model system reliability. This method builds on the Piecewise NHPP model by introducing the concept of change points, which represent the time when the underlying failure behavior changes. Detection of change points involves statistical techniques that analyze failure data to automatically identify when the behavior changes, allowing for a more precise segmentation of the model into different failure distributions.

### Example

To run a Piecewise NHPP with Change Point Detection with ReliaGrowR, first, set up some cumulative time and failure data:

```{r}
times <- c(25, 55, 97, 146, 201, 268, 341, 423, 513, 609, 710, 820, 940, 1072, 1217)
failures <- c(1, 1, 2, 4, 4, 1, 1, 2, 1, 4, 1, 1, 3, 3, 4)
```

Then run the rga and plot the results:

```{r}
result <- rga(times, failures, model_type = "Piecewise NHPP")
plot_rga(result)
```

## Accelerated Life Testing

**Accelerated Life Testing (ALT)** is a method to estimate the life of a product or system under normal operating conditions by subjecting it to accelerated stress levels in a controlled environment. The goal is to identify potential failure modes and predict how long the product will last under normal conditions without having to wait for failures to occur naturally over a longer time period.

### Accelerated Life Testing: Stress vs Life

Take a product that is designed to last 10,000 hours. Waiting for it to fail naturally could take months. With ALT, you might increase the operating temperature of the product, causing it to fail faster. 

The following plot shows the relationship between the acceleration or stress factor and the expected life of the product: as the acceleration factor increases, the expected life tends to decrease. Adjust the acceleration factor input to see this impact in action. 

```{r echo=FALSE}
sliderInput("af", label = h3("Acceleration Factor:"), value = 1, min = 1, max = 10, step = 0.1, width = '50%')
plotOutput("altPlot")
```

```{r, context="server"}
output$altPlot <- renderPlot({
  y <- seq(1, 10, by = .1)
  x <- 3678.79 * exp(1 / y)
  life <- 3678.79 * exp(1 / input$af)
  plot(x, y, type = "l", col = "blue", xlab = "Life", ylab = "Stress Factor", main = "Stress Factor vs Life")
  abline(v = life, col = "red", lty = 2)
  legend("topright", legend = c("Expected Life = ", ceiling(life)))
})
```

## The Arrhenius Model

The **Arrhenius Model** is widely used in ALT to model the effect of temperature on the rate of chemical reactions, which directly relates to the degradation or failure rate of materials.

The Acceleration Factor (AF) for the Arrhenius model is given by:

$$ AF = \exp\left(\frac{E_a}{k} \left( \frac{1}{T_{\text{use}}} - \frac{1}{T_{\text{stress}}} \right) \right) $$

Where:

- \( E_a \) = Activation energy in eV
- \( k \) = Boltzmann constant (\( 8.617 \times 10^{-5} \, \text{eV/K} \))
- \( T_{\text{use}} \) = Temperature at normal use condition in Kelvin
- \( T_{\text{stress}} \) = Temperature at accelerated stress condition in Kelvin

By fitting this model to test data, you can estimate the AF and predict product life under normal operating conditions.

## WeibullR.ALT

For this tutorial, we will use **WeibullR.ALT**, an extension of the WeibullR package for Accelerated Life Testing.

First, check if WeibullR.ALT is installed in R and install if not.

```{r echo=TRUE, message=FALSE, warning=FALSE}
if(!(require(WeibullR.ALT))){install.packages('WeibullR.ALT')}
```

For this example, we will use the data set `MeekerData` that comes prepackaged with `WeibullR.ALT`. First, load the data set and view the first few samples.

```{r}
data <- MeekerData("table10")
head(data)
```

The data set contains columns for the time to failure **time**, event type **event** (1 = failure, 0 = suspension), quantity of samples **qty** and temperature as the acceleration factor **TempC**. Testing was conducted at 4 different temperature levels, 10, 40, 60, and 80 degrees C. 

Next, use the `alt.data` function to create a separate object for the data at each temperature level. 

```{r}
data10C <- alt.data(data[which(data$TempC==10), 1:3], 
                    stress = data[which(data$TempC==10), 4][1])
data40C <- alt.data(data[which(data$TempC==40), 1:3], 
                    stress = data[which(data$TempC==40), 4][1])
data60C <- alt.data(data[which(data$TempC==60), 1:3], 
                    stress = data[which(data$TempC==60), 4][1])
data80C <- alt.data(data[which(data$TempC==80), 1:3], 
                    stress = data[which(data$TempC==80), 4][1])
```

Then, use the `alt.make` function to fit an Arrhenius-Weibull model to the data.

```{r}
arrobj <- alt.make(list(data10C, data40C, data60C, data80C), dist = "weibull",
                alt.model = "arrhenius", method.fit = "mle")
```

Finally, calculate the AF at 80 degrees C using the Arrhenius Model. Assume that the Activation Energy is 0.7 eV.

```{r, results='asis'}
E_a <- 0.7
k <- 8.617E-5
T_use <- 10 + 273
T_stress <- 80 + 273
AF <- exp(E_a / k * (1 / T_use - 1 / T_stress))
cat(paste0("The Acceleration Factor at 80 degrees C is ", round(AF, 2)))
```

### Exercises

```{r arrhns1, echo=FALSE}
question("What is the approximate unreliability for a duration of 3000 and an operating stress of 60?",
  answer("90%"),
  answer("10%"),
  answer("70%"),
  answer("30%", correct = TRUE),
  random_answer_order = TRUE
)
```

```{r arrhns2, echo=FALSE}
question("What is the approximate reliability for a duration of 1000 and an operating stress of 80?",
  answer("70%"),
  answer("30%"),
  answer("40%"),
  answer("60%", correct = TRUE),
  random_answer_order = TRUE
)
```

## The Power Law Model

The **Power Law Model** is another commonly used model in ALT, particularly for mechanical systems and components. This model describes how a product's life is affected by a certain stress factor (such as temperature, voltage, or pressure). Unlike the Arrhenius model, which is primarily used for temperature-based stresses, the Power Law Model can be applied to a wide variety of stress factors.

The AF for the Power Law model is given by:

$$ AF = \left( \frac{S_{\text{stress}}}{S_{\text{use}}} \right)^{n} $$

Where:

- \( S_{\text{stress}} \) = Stress level during accelerated testing (e.g., voltage, load)
- \( S_{\text{use}} \) = Stress level during normal use conditions
- \( n \) = Stress exponent, indicating the sensitivity of the failure mechanism to the applied stress

For this example, we will use the data set `NelsonData` that comes prepackaged with `WeibullR.ALT`. First, load the data set and view the class of data.

```{r results='asis'}
data <- NelsonData("table3.1")
cat(class(data))
```

The data is a list of time to failure data at 7 different voltage stress levels, 26KV, 28kV, 30kV, 32kV, 34kV, 36kV, and 38kV.

Next, use the `alt.data` function to create a separate object for the data at each stess level. 

```{r}
data26 <- alt.data(data$kV26, stress = 26)
data28 <- alt.data(data$kV28, stress = 28)
data30 <- alt.data(data$kV30, stress = 30)
data32 <- alt.data(data$kV32, stress = 32)
data34 <- alt.data(data$kV34, stress = 34)
data36 <- alt.data(data$kV36, stress = 36)
data38 <- alt.data(data$kV38, stress = 38)
```

Then, use the `alt.make` function to fit an Power-Weibull model to the data.

```{r}
pwrobj <- alt.make(list(data26, data28, data30, data32, data34, data36, data38), 
                   dist = "weibull", alt.model = "power")
```

### Exercises

```{r pwr1, echo=FALSE}
question("Which subset of data does not seem to follow the Power Law model?",
  answer("The data with a stress level of 38kV"),
  answer("The data with a stress level of 28kV"),
  answer("The data with a stress level of 36kV"),
  answer("The data with a stress level of 26kV", correct = TRUE),
  random_answer_order = TRUE
)
```

```{r pwr2, echo=FALSE}
question("How would you interpret a plot to the right of another plot (e.g. set 2 in green vs set 3 in red)?",
  answer("The plot to the right has a lower unreliability than the one to the left", correct = TRUE),
  answer("The plot to the right a higher unreliability than the one to the left"),
  answer("The plot to the right has a lower reliability than the one to the left"),
  answer("The plot to the right has a higher reliability than the one to the left", correct = TRUE),
  random_answer_order = TRUE
)
```

## Multiple Stress Factors

The **Eyring Model** is an extension of the Arrhenius model and is commonly used in ALT to account for the combined effects of multiple stressors such as temperature, voltage, and humidity on the failure mechanisms of a product. This model is particularly useful when the product is exposed to environmental factors beyond just temperature.

The general form of the Eyring Model is:

\[\text{Life}(S) = A \cdot \exp\left(\frac{B}{S} + C \cdot S\right)\]

Where:

* \( \text{Life}(S) \) is the expected life under stress level \( S \),
* \( A \) is a constant,
* \( B \) and \( C \) are model parameters,
* \( S \) is the applied stress factor (e.g., temperature, voltage, humidity).

## Parallellization

A common assumption in ALT is that the subset of data at each stress level have the same slope or shape factor. If each model subset has the same slope, then all models appear as parallel plots. 

To fit a parallel model to the Meeker data, use the `alt.parallel` function.
```{r}
prlobj <- alt.parallel(arrobj)
```

Similarly, a parallel model can be fitted to the Nelson data.
```{r}
expobj <- alt.parallel(pwrobj)
```

## Relationship Plots

Another way to view the relationship between the stress factor and the life of a product is with a relationship plot. A relationship plot is a log-linear plot with the stress factor on the x-axis and the log of the unreliability (failure probability) on the y-axis.

A linear model is fit to the data, which shows patterns indicating how life changes as stress increases. 

* A **negative slope** (decreasing life with increasing stress) indicates that higher stress leads to shorter life. 
* A **positive slope** (increasing life with increasing stress) is less common but indicates that life improves with higher stress.
* A **flat slope** indicates that higher stress does not have an impact on life. 

To fit a relationship plot to the Meeker data, first fit a linear model to the data and then plot the results.
```{r}
lnrobj <- alt.fit(prlobj)
plot(lnrobj, suppress.dev.new = TRUE)
```

Similarly, a relationship plot can be fitted to the Nelson data.

```{r}
lnrobj2 <- alt.fit(expobj)
plot(lnrobj2, suppress.dev.new = TRUE)
```

### Exercises

```{r rel1, echo=FALSE}
question("What does the slope of the power-Weibull relationship plot indicate?",
  answer("Higher stress does not seem to have an impact on life"),
  answer("Higher stress tends to lead to longer life"),
  answer("Higher stress tends to lead to shorter life", correct = TRUE),
  random_answer_order = TRUE
)
```

```{r rel2, echo=FALSE}
question("On the power-Weibull relationship plot, what is the approximate P90 time to failure at a stress level of 30?",
  answer("100", correct = TRUE),
  answer("10"),
  answer("1000"),
  answer("1"),
  random_answer_order = TRUE
)
```

## Summary

In this tutorial, we've introduced the concepts of Reliability Growth Analysis and Accelerated Life Testing. You have learned to calculate and visualize reliability using different types of models.

### More resources

[WeibullR.plotly](https://paulgovan.github.io/WeibullR.plotly/) is a package for building interactive Weibull models. 

[WeibullR.shiny](https://paulgovan.github.io/WeibullR.shiny/) is a web application for life data analysis.

### References

- Abernethy, R.B. (2004). The New Weibull Handbook. Fifth Edition.

- Aden-Buie G, Schloerke B, Allaire J, Rossell Hayes A (2023). learnr: Interactive Tutorials for R. https://rstudio.github.io/learnr/, https://github.com/rstudio/learnr.

- Crow, Larry H. (1975). Reliability Analysis for Complex, Repairable Systems. Army Material Systems Analysis Activity, 40. <https://apps.dtic.mil/sti/citations/ADA020296>

- Fasola S, Muggeo VMR, Kuchenhoff K. (2018). A heuristic, iterative algorithm for change-point detection in abrupt change models. Computational Statistics, 33, 997-1015.

- Govan P (2024). ReliaGrowR: Reliability Growth Analysis. R package version 0.1,
<https://cran.r-project.org/package=ReliaGrowR>.

- Guo, H., Mettas, A., Sarakakis, G. and Niu, P., (2010). Piecewise NHPP models with maximum likelihood estimation for repairable systems. 2010 Proceedings - Annual Reliability and Maintainability Symposium (RAMS), 1-7. <https://doi.org/10.1109/RAMS.2010.5448029>.
  
- Duane J. T. (1964). Learning Curve Approach to Reliability Monitoring. IEEE Transactions on Aerospace, 2, 563-566, <https://doi.org/10.1109/TA.1964.4319640>.

- Meeker W. Q., Escobar L. A., (1998). Statistical Methods for Reliability Data, Wiley-Interscience, New York.

- Nelson W, (1990). Accelerated Testing, Wiley-Interscience, New York.

- Silkworth, David. (2020). WeibullR: An R Package for Weibull Analysis for Reliability Engineers. 43-53. <https://doi.org/10.35566/isdsa2019c3>.

- Silkworth D, Symynck J (2022). WeibullR: Weibull Analysis for Reliability Engineering. R package version 1.2.1, <https://CRAN.R-project.org/package=WeibullR>.

- Silkworth D (2022). WeibullR.ALT: Accelerated Life Testing Using 'WeibullR'. R package version 0.7.2, <https://CRAN.R-project.org/package=WeibullR.ALT>.
