Reliability Analysis for Degradation-Shock Processes with State-Varying Degradation Patterns Using Approximate Bayesian Computation (ABC) for Parameter Estimation

Abstract

The degradation of products is an integral part of their life-cycle, often following predictable trajectories. However, sudden, unexpected events, termed ’shocks’, can substantially alter these degradation paths. Shocks can significantly influence the pace of degradation, leading to accelerated system failure. Moreover, they may initiate changes in degradation patterns, transitioning from linear to non-linear or random trajectories. To address this challenge, we present a novel multi-state reliability model for competing failure processes that account for degradation-shock dependencies by considering the state-varying degradation pattern. The degradation process is divided into s-states, with each state treated according to its pattern based on the time-transform Wiener process. The reliability function is derived based on soft failure caused by continuous degradation involving the s-states, the sudden increase in degradation caused by random shocks, and hard failure due to some shock processes. Additionally, we performed a sensitivity analysis to determine which parameters have the most significant impact on product reliability. Due to the complexity of the likelihood function, we adopted the ABC method to estimate the model parameters. A simulation study and a practical application with micro-electro-mechanical systems (MEMS) degradation results are used to demonstrate the efficiency and effectiveness of the proposed approach.

1. Introduction

Over the past few years, there has been a significant increase in the use of degradation data to assess the reliability of products. These data provide essential insights into a product’s health and can be obtained through non-destructive and non-invasive methods, primarily by monitoring sensors [1,2]. This approach is not only cost-effective but also more convenient than destructive testing. Degradation analysis can reveal underlying failure mechanisms and assist in forecasting a product’s future performance. Consequently, degradation analysis has gained widespread acceptance as a means of studying product reliability such us the failure degradation of lithium battery [3,4,5], train wheel wear degradation [6], and so forth. There are various approaches to modeling degradation data, including the general path method and stochastic process models [7]. However, stochastic process models have gained popularity in degradation modeling due to their strong mathematical properties and their ability to accurately represent the stochastic nature of degradation data [8,9]. Several widely used models in this domain include the Wiener process [10,11], Gamma process [12,13], and inverse Gaussian process [14,15], among others.

In many real-world engineering scenarios, the reliability of a product is affected by various factors such as environmental conditions, wear and tear, and other forms of degradation. A product may exist in different states of degradation throughout its life cycle, and each state may have a different level of reliability and performance. Therefore, multi-state reliability modeling is an effective approach to analyzing complex products that exhibit multiple degradation states. Random shocks, also known as shock loads, can significantly contribute to product failures and impact the degradation process. These shocks can cause stress concentrations and damage to components, leading to accelerated degradation and ultimately resulting in product failure [16]. It is, therefore, crucial to consider the correlation between the shock process and the degradation process when developing effective reliability models for complex engineering products [17]. Typically, a product can fail due to two failure modes: soft failure from degradation or hard failure from random shocks. Regardless of the failure mode, it can result in the product’s failure [18]. Figure 1 illustrates the degradation-shock competing failure process of the product.

The development of models for systems experiencing competing failure processes has long captured the attention of researchers. The foundational work of Lemoine and Wenocur [19] marked the early exploration of this complex issue. Over the years, numerous researchers have advanced this field, proposing various models to address multiple failure processes. For instance, Peng et al. [20] developed reliability models for complex systems that experience multiple dependent competing failure processes, taking into consideration the failures caused by degradation and shock loads. Jiang et al. [21] developed reliability models for complex systems experiencing multiple dependent competing failure processes, accounting for failures caused by degradation and shock loads, with a dynamic failure threshold that decreases over time due to shocks. Song et al. [22] developed a reliability model for complex multi-component systems where each component experiences multiple dependent failure processes due to degradation and shared shock loads. Fan et al. [23] developed a reliability model for dependent competing failure processes with degradation-shock dependence, classifying shocks into damage, fatal, and safety zones. Che et al. [24] developed a reliability model for systems with mutually dependent degradation and shock processes. They introduced a facilitation model, a Markov point process, to capture the shock process, which cannot be described by traditional Poisson models. Cao and Dong [25] analyzed the competing and dependent failure processes for multi-state systems suffering from four typical random shocks. Liang et al. [26] developed a reliability model for a system subject to both multi-state and continuous degradation processes. In their work, the degradation process was modeled by a time-homogeneous semi-Markov process on a finite state space and a time-homogeneous stochastic process on a continuous state space with independent increments. Feng et al. [27] proposed a degradation-shock-dependent competing failure process model considering a time-shifting sudden failure threshold for predicting the RUL of drill bits. Chang et al. [28] proposed a generalized reliability model for systems with arbitrary structures experiencing both degradation and shock processes. Shao et al. [29] proposed a hybrid RUL prediction method for a subsea hydraulic control system that considers degradation-shock dependency, introducing a dependency factor to model the interaction between shock intensity and degradation stages. Lyu et al. [30] developed a reliability model for systems subject to dependent competing failure processes, incorporating the effects of random cycle times, where degradation rates and hard failure thresholds vary with the number of cycles. Qiu et al. [31] derived the recursive equation for systems reliability and availability considering the impact of shock in the degradation process and the catastrophic failure. Jin and Zhang [32] explored cascading failures in circuit systems by modeling the system as an impedance network. Their research emphasized the coupling effects of continuous degradation and random shocks, showing how these interactions can lead to asymmetric failure propagation and unpredictable failure times. Gan and Tang [33] developed a reliability model that accounts for the degradation shock process, considering the dynamic nature of failure thresholds. They used a linear Wiener process to describe performance degradation and introduced a multi-performance degradation model based on copula functions. Xu et al. [34] proposed a multi-stage model that accounts for the dependencies between degradation and shock processes across different damage levels.

Notably, most of these mentioned studies considered that, during the soft failure process, the arrival of shock may result in accelerating the degradation process. However, it is possible that the arrival of a shock may not only accelerate the degradation process but also change the degradation pattern. For example, the shock may cause changes in the distribution of stresses and strains within the material, leading to a shift in the location and nature of degradation mechanisms [35]. Additionally, the shock may induce new degradation mechanisms that were not present before, leading to a change in the degradation pattern. In some cases, a shock may cause a sudden increase in the degradation rate, followed by a gradual return to the pre-shock degradation pattern [36]. In other cases, the degradation pattern may permanently shift from linear to non-linear, or from one type of non-linear behavior to another [37], for example, a mechanical component such as metal gear that is subjected to cyclic loading. Under normal operating conditions, the gear undergoes linear fatigue degradation, meaning that the number of cycles until failure decreases linearly over time. However, if the gear is subjected to an unexpected shock, such as a sudden increase in the load or a change in the loading direction, the degradation pattern may change to include non-linear fatigue degradation. In this case, the number of cycles until failure may decrease at an accelerated rate, exhibiting a more rapid degradation pattern than the linear fatigue model would predict [38]. Similarly, let us consider an electrical component, such as a capacitor. Under normal operating conditions, the capacitor undergoes linear degradation due to electrolyte drying, meaning that the capacitance decreases linearly over time. However, if the capacitor is subjected to an unexpected shock, such as a sudden increase in the operating temperature or a change in the applied voltage, the degradation pattern may change to include non-linear degradation due to dielectric breakdown. In this case, the capacitance may decrease more rapidly than predicted by the linear model, exhibiting a more complex degradation pattern [39,40].

It is less common for a shock arrival to change the degradation pattern from non-linear to linear, as shocks typically introduce additional stress and damage to a component, which can lead to more complex degradation mechanisms. However, one scenario where this could potentially occur is if a component is undergoing a complex degradation process involving multiple mechanisms, and the shock arrival causes some of the mechanisms to become temporarily inactive. For example, consider a metal component that is undergoing wear and corrosion degradation due to exposure to a harsh environment. Under normal operating conditions, the degradation may occur in a non-linear fashion due to the interplay of the mechanical and chemical degradation mechanisms. However, if the component experiences a sudden shock or impact, it could cause the corrosion mechanism to become temporarily inactive due to the mechanical damage caused by the shock. This could lead to a more linear degradation pattern dominated by wear, at least until the corrosion mechanism becomes reactivated. Therefore, it’s essential to consider these dynamic impact of shocks when developing models for multi-state degradation-shock failure processes. Figure 2 depicts common degradation-shock processes that often appear in practice.

In terms of parameter estimation, the ABC method is selected due to the complexity of our likelihood function. This algorithm approximates the posterior distribution $\pi \left(\theta \mid X_{\mathrm{obs}}\right)$ by first generating parameter samples from the joint distribution of the parameter vector $\theta =\left\{{\theta }_1,\dots ,{\theta }_n\right\}$, then simulating pseudo-data based on these sampled parameters from the model $\mathcal{M}\left(\theta \right)$, and finally accepting those parameter samples that produce pseudo-data close to the observed data, satisfying $\rho (X_{\mathrm{obs}},X)\le \epsilon$. When working with high-dimensional data, it is more efficient to define the distance measure based on summary statistics of both the observed and simulated data: $\rho (S\left(X_{\mathrm{obs}}\right),S\left(X\right))\le \epsilon$, where $\rho (.)$ is the distance function, $\epsilon$ is the chosen tolerance level, and $S(.)$ refers to the summary statistics. A smaller $\epsilon$ results in a better approximation. For further details, refer to [41,42].

The ABC method is highly dependent on the selection of a distance measure to evaluate the difference between observed and simulated data. Various distance measures, such as Manhattan [43], Euclidean [44], Canberra [45], and Chebyshev [46], have been commonly used, but the optimal choice depends on the nature of the data. In this study, the Euclidean distance was applied due to its simplicity and ease of interpretation. This distance is often used in ABC to measure differences between the summary statistics of simulated and observed data rather than comparing entire datasets, especially when dealing with complex, high-dimensional data. In our work, the mean was chosen as the primary summary statistic, justified by the Wiener process that underlies our degradation model. Since this process involves normally distributed increments, using the mean effectively captures the central trend of the data, making it an optimal and computationally straightforward choice. To achieve convergence in the ABC algorithm, careful attention was given to both the summary statistics and the calibration of the tolerance threshold $\epsilon$. This threshold was chosen to strike a balance between the accuracy of parameter estimation and computational efficiency. If the threshold is too low, the computational cost rises without significant improvement in accuracy, particularly when the summary statistics do not fully represent the variability in the data. The ideal threshold depends on factors like data variability and simulation acceptance rates, ensuring that the algorithm converges effectively to the true parameter values without overfitting.

In recent years, various sampling algorithms like Markov-chain Monte Carlo (MCMC), population Monte Carlo, and sequential Monte Carlo have been combined with ABC to enhance its accuracy and efficiency for parameter estimation. These extended methods include ABC-MCMC [42,47], ABC-HMC [48], ABC-PMC [49], and ABC-Gibbs [50,51], with ABC-Gibbs being especially useful for high-dimensional models.

ABC has gained popularity across fields such as population genetics and engineering. For instance, Rodrigues et al. used ABC-Gibbs to estimate parameters in a high-dimensional state-space model for Airbnb rental prices, and Muhammad et al. applied it to a complex degradation model. Based on these successes, we used ABC-Gibbs in our work for parameter estimation. More details on ABC-Gibbs can be found in [50,51]. The main contributions of this paper can be summarized as follows:

• We propose a multi-state reliability model that accounts for the degradation-shock dependency and considers state-varying degradation patterns using a time-transform Wiener process. Our approach is based on dividing the degradation process into s-states and treating each state according to its pattern based on the time-transform Wiener process.
• We derive the reliability function considering soft failure caused by continuous degradation across multiple states, sudden increases in degradation due to random shocks, and hard failure from certain shock processes.
• We adopt approximate Bayesian computation (ABC) for parameter estimation in complex reliability models, overcoming the limitations posed by highly complex likelihood functions.
• We perform a comprehensive sensitivity analysis to identify the most influential parameters affecting system reliability, providing critical insights for system design and maintenance.

The remainder of the paper is organized as follows: Section 2 outlines the reliability analysis framework for systems affected by random shocks and their impact on degradation trajectories. Section 3 presents a simulation study and a practical application of the model to micro-electro-mechanical systems (MEMS) degradation data, illustrating the effectiveness and efficiency of our approach. Finally, Section 4 provides the conclusion of the study.

2. Reliability Modeling of a System Subjected to a Degradation-Shock Process with Multiple States

This section focuses on the reliability analysis of systems and components considering the influence of random shocks and the resulting changes in the degradation pattern. Shocks can accelerate the degradation process and lead to shifts from linear to non-linear or random degradation trajectories. The objective is to develop precise expressions for the reliability function, accounting for extreme shocks and the multi-phased nature of degradation.

2.1. Shock Modeling Based on Extreme Shock

Under the extreme shock process, the probability that the system has not failed by time t can be calculated as the probability that the degradation level remains below the hard threshold value $X_w$ at that time. Specifically, T denotes the random time of hard failure, occurring when cumulative shocks exceed the threshold $X_w$. The expression $P(t\le T)$ thus calculates the probability of survival, indicating that degradation due to shocks up to time t remains under the critical level. Assuming that the magnitudes of shocks $W_i$ are independent and identically distributed random variables following a normal distribution with mean ${\mu }_w$ and variance ${\sigma }_w^2$, this probability can be expressed as:

$$P(t\le T)=P\left(W_1<X_w,W_2<X_w,\dots ,W_{N\left(t\right)}<X_w\right).$$

(1)

The random variable $N\left(t\right)$ represents the number of shocks that occur in the time interval $[0,t]$, and is modeled as a Poisson process with rate $\lambda$. The probability of $N\left(t\right)$ shocks occurring in the time interval $[0,t]$ is given by the Poisson distribution:

$$P(N\left(t\right)=s)=e^{-\lambda t}{\displaystyle \frac{{\left(\lambda t\right)}^s}{s!}},s=0,1,2,\dots .$$

(2)

Given that $N\left(t\right)$ shocks have occurred in the time interval $[0,t]$, the probability of all shocks having magnitudes less than the threshold value $X_w$ is given by:

$$P(W_1<X_w,W_2<X_w,\dots ,W_{N\left(t\right)}<X_w)={\left[\Phi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\right]}^{N\left(t\right)}.$$

(3)

where $W_i\ge 0,{\mu }_w>0,{\sigma }_w$ is assumed to be nearly zero to minimize the occurrence of negative values, and $\Phi (·)$ is the cumulative distribution function of the standard normal distribution.

2.2. Degradation Modeling under Soft Failure and Pattern Variation

Degradation processes are ubiquitous in numerous systems, leading to a gradual decline in their performance over time. Factors such as wear, corrosion, fatigue, and aging, among others, trigger these degradation processes, ultimately impacting the reliability, functionality, and safety of the system. Therefore, accurately estimating the lifetime of a system based on degradation data is essential for effective maintenance, decision making, and resource allocation [10]. The precision of this life estimate significantly relies on the degradation model’s ability to capture the complexities of the system’s failure degradation process. Additionally, random shocks not only accelerate degradation but also alter its pattern, potentially resulting in soft failure once the soft failure threshold is exceeded. To incorporate this scenario into the degradation model, the degradation process is divided into s distinct states, each characterized by its unique degradation pattern. In addition to that, each state is modeled by the time-transform Wiener process. The “time-transform Wiener process” refers to a Wiener process where time is transformed using a function ${\Lambda }_s\left(t\right)=t^{b_s}$, allowing for flexible modeling of different degradation rates across states. More detail regarding the time transform Wiener process can be found in [10,52].

Suppose that the degradation process comprises s states, with each state being initiated due to shock arrival, lasting within the time interval $[t_s,t_{s+1}]$. The degradation process for each state can then be expressed as:

$$\begin{array}{ccc}\hfill X_0\left(t\right)& =x_{t_0}+{\mu }_0{\Lambda }_0\left(t\right)+{\sigma }_{0}B\left({\Lambda }_0\left(t\right)\right)\hfill & \hfill t_0\le t\le t_1\\ \hfill X_1\left(t\right)& =x_{t_1}+{\mu }_1{\Lambda }_1\left(t\right)+{\sigma }_{1}B\left({\Lambda }_1\left(t\right)\right)\hfill & \hfill t_1\le t\le t_2\\ \hfill X_2\left(t\right)& =x_{t_2}+{\mu }_2{\Lambda }_2\left(t\right)+{\sigma }_{2}B\left({\Lambda }_2\left(t\right)\right)\hfill & \hfill t_2\le t\le t_3\\ \hfill & ⋮\hfill \\ \hfill X_s\left(t\right)& =x_{t_s}+{\mu }_s{\Lambda }_s\left(t\right)+{\sigma }_{s}B\left({\Lambda }_s\left(t\right)\right)\hfill & \hfill t_s\le t\le t_{s+1}\end{array}$$

(4)

where $x_{t_0}$ is assumed to be zero for simplicity and $X_s\left(t\right)$ represents the degradation process of s-th state, which evolves over time t between $t_s$ and $t_{s+1}$. The initial value of $X_s$ at time $t_s$ is given by $x_{t_s}$. The time interval $[t_s,t_{s+1}]$ represents the time span during which the s-th state evolves according to this equation. Note that the upper limit of this interval, $t_{s+1}$, is assumed to be the same as the lower limit of the next interval. ${\mu }_s$ and ${\sigma }_s$ are the drift and diffusion coefficients of s-th state and $B\left({\Lambda }_s\left(t\right)\right)$ is a standard Brownian motion process. ${\Lambda }_s\left(t\right)$ represents a transformed time scale for the s-th state, defined as ${\Lambda }_s\left(t\right)=t^{b_s}$, where $b_s$ is a parameter determining the nature of the time transformation. If $b_s=1$, the transformation is linear, indicating that the degradation progresses at a constant rate over time. However, if $b_s\ne 1$ then the process is non-linear; this flexibility allows the model to capture the varying dynamics of degradation under different conditions, especially in response to shock arrivals, where each state’s degradation pattern can change significantly.

The shocks cause abrupt jumps of sizes $h_1,h_2,\dots ,h_{N\left(t\right)}$, which results in increasing the degradation process. Let $H\left(t\right)$ be the cumulative length of the shock during the time interval $(0,t]$, then $H\left(t\right)$ can be expressed as:

$$H\left(t\right)=\left\{\begin{array}{cc}\sum _{j=1}^{N\left(t\right)}{h}_j,\hfill & \mathrm{if}N\left(t\right)>0\hfill \\ 0,\hfill & \mathrm{if}N\left(t\right)=0.\hfill \end{array}\right.$$

(5)

Therefore, the total length of the degradation path can be expressed as

$$X\left(t\right)=\sum _{i=0}^s{X}_i\left(t\right)+\sum _{j=1}^{N\left(t\right)}{h}_j$$

(6)

where $N\left(t\right)$ indicates the number of shocks. Based on the properties of the Wiener process, $X_i\left(t\right)$ follows normal distribution with mean $x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))$ and variance ${\sigma }_i{\Lambda }_i\left(t\right)$, i.e $X_i\left(t\right)\sim N(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right)),{\sigma }_i{\Lambda }_i\left(t\right))$ for $1\le i\le s$. If $h_i$ is independent and identically distributed (i.i.d), and follows normal distribution with mean ${\mu }_h$ and variance ${\sigma }_h^2$, i.e $h_i\sim N({\mu }_h,{\sigma }_h^2)$ for $1\le i\le N\left(t\right)$. Then, $X\left(t\right)$ follows normal distribution, $X\left(t\right)\sim N\left({\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+{\mu }_{h}N\left(t\right),{\sum }_{i=1}^s{\sigma }_i^2{\Lambda }_i\left(t\right)+N\left(t\right){\sigma }_h^2\right)$.

If we assumed that there exist uncertain errors during the process of data collection, then the length of the degradation path could be modified as

$$\begin{array}{cc}\hfill & X_{e_0}\left(t\right)=X_0\left(t\right)+{ϵ}_0\hfill \\ \hfill & X_{e_1}\left(t\right)=X_1\left(t\right)+{ϵ}_1\hfill \\ \hfill & X_{e_s}\left(t\right)=X_s\left(t\right)+{ϵ}_s\hfill \\ \hfill & ⋮\hfill \\ \hfill & X_e\left(t\right)=\sum _{i=0}^s{X}_{e_i}\left(t\right)\hfill \end{array}$$

(7)

where ${ϵ}_0,{ϵ}_1,\dots ,{ϵ}_s$, represents independent and identical measurement errors and are assumed to follow a normal distribution with mean zero and standard deviation ${\sigma }_{ϵ}$, i.e., ${ϵ}_i\sim N(0,{\sigma }_{ϵ}^2)$. Hence, $X_e\left(t\right)\sim N\left({\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+{\mu }_{h}N\left(t\right),{\sum }_{i=1}^s{\sigma }_i^2{\Lambda }_i\left(t\right)+N\left(t\right){\sigma }_h^2+s{\sigma }_{ϵ}^2\right)$.

2.3. Reliability Modeling

The product fails when either of the two failure modes occurs—that is, soft failure or hard failure. To model the reliability under the influence of shock and the degradation process, we assume that the system experiences a series of shocks at random intervals of time $[t_s,t_{s+1}]$. We also assume that the degradation state changes from state s to state $s+1$ and that, at each state, the degradation pattern is considered changed. $X_f$ denotes the failure threshold of the degradation process, representing the critical level at which the cumulative degradation is considered sufficient to cause system failure due to soft failure. The system is assumed to have failed when the degradation process reaches or exceeds this threshold. Hence, the probability of the system at state i can be expressed as:

$$\begin{array}{cc}\hfill P_0& =P_r\left(N\left(t\right)=0,X_e\left(t\right)<X_f\right)=P_r(N\left(t\right)=0,X_0\left(t\right)+{ϵ}_0<X_t)\hfill \\ & =P_r(N\left(t\right)=0)P_r(X_0\left(t\right)+{ϵ}_0<X_f),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill P_1& =P_r\left(N\left(t\right)=1,X_e\left(t\right)<X_f,W_1<X_w\right)\hfill \\ \hfill & =P_r\left(N\left(t\right)=1,X_0\left(t\right)+X_1\left(t\right)+h_1+{ϵ}_0+{ϵ}_1<X_f,W_1<X_w\right)\hfill \\ \hfill & =P_r(N\left(t\right)=1)P_r\left(X_0\left(t\right)+X_1\left(t\right)+h_1+{ϵ}_0+{ϵ}_1<X_f\right)P_r\left(W_1<X_w\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill P_2=& P_r\left(N\left(t\right)=2,X_e\left(t\right)<X_f,W_1<X_w,W_2<X_w\right)\hfill \\ \hfill =& P_r\left(N\left(t\right)=2,X_0\left(t\right)+X_1\left(t\right)+X_2\left(t\right)+h_1+h_2+{ϵ}_0+{ϵ}_1+{ϵ}_2<X_f,W_1<X_w,W_2<X_w\right)\hfill \\ \hfill =& P_r(N\left(t\right)=2)P_r\left(X_0\left(t\right)+X_1\left(t\right)+X_2\left(t\right)+h_1+h_2+{ϵ}_0+{ϵ}_1+{ϵ}_2<X_f\right)P_r\hfill \\ \hfill & \left(W_1<X_w\right)P_r\left(W_2<X_w\right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill P_{\mathrm{s}}=& P_r\left(N\left(t\right)=\mathrm{s},X_e\left(t\right)<X_f,W_1<X_w,W_2<X_w,\cdots ,W_s<X_w\right)\hfill \\ \hfill & =P_r(N\left(t\right)=\mathrm{s})\hfill \\ \hfill & P_r\left(X_1\left(t\right)+X_2\left(t\right)+X_3\left(t\right)+\cdots +X_s\left(t\right)+h_1+h_2+\cdots +h_s+{ϵ}_0+{ϵ}_1+\cdots +{ϵ}_i<X_f\right)\hfill \\ \hfill & P_r\left(W_1<X_w\right)P_r\left(W_2<X_w\right)\cdots P_r\left(W_n<X_w\right)\hfill \\ \hfill =& P_r(N\left(t\right)=s)\times P_r\left(\sum _{i=0}^s{X}_i\left(t\right)+\sum _{i=1}^s{h}_i+\sum _{i=0}^s{ϵ}_i<X_f\right)\times \prod _{i=1}^s{P}_r(W_i<X_w)\hfill \\ \hfill =& P_r(N\left(t\right)=s)\times P_r\left(\sum _{i=0}^s{X}_i\left(t\right)+\sum _{i=1}^s{h}_i+\sum _{i=0}^s{ϵ}_i<X_f\right){\left[P_r\left(W_1<X_w\right)\right]}^s\hfill \end{array}$$

(11)

Therefore, the reliability of the system under extreme shock can be expressed as

$$\begin{array}{cc}\hfill R_e\left(t\right)=\sum _{i=0}^s{P}_i=& \Phi \left({\displaystyle \frac{X_f-{\mu }_0{\Lambda }_0\left(t\right)}{\sqrt{{\sigma }_0^2{\Lambda }_0\left(t\right)+{\sigma }_{ϵ}^2}}}\right)exp(-\lambda t)+\sum _{i=1}^s{\left[\Phi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\right]}^i\hfill \\ \hfill & \Phi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+i{\sigma }_{ϵ}^2}}}\right){\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}.\hfill \end{array}$$

(12)

2.4. Reliability Sensitivity Analysis

Sensitivity analysis is a valuable tool for assessing the influence of various parameters and assumptions on the reliability model. By conducting sensitivity analysis, we can gain insights into which factors have the most significant impact on the reliability function and better understand the robustness of the model. In this section, we will perform a sensitivity analysis for the reliability model based on extreme shocks and multi-phased degradation process.

2.4.1. Sensitivity of Reliability with Respect to Shock Parameters

The sensitivity of the reliability with respect to shock parameters can be determined by differentiating the reliability function with respect to each parameter. Considering the shock parameters $\lambda$, ${\mu }_w$, and ${\sigma }_w$, the sensitivities of reliability $R_e\left(t\right)$ can be computed as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial \lambda }}& =-t\Phi \left({\displaystyle \frac{X_f-{\mu }_0{\Lambda }_0\left(t\right)}{\sqrt{{\sigma }_0^2{\Lambda }_0\left(t\right)+{\sigma }_{\epsilon }^2}}}\right)exp(-\lambda t)\hfill \\ \hfill & -\sum _{i=1}^s{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\Phi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\mu }_w}}& =-\sum _{i=1}^{s}i{\left[\Phi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\right]}^{i-1}{\displaystyle \frac{1}{{\sigma }_w}}\varphi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\sigma }_w}}& =-\sum _{i=1}^{s}i{\left[\Phi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\right]}^{i-1}\left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w^2}}\right)\varphi \left({\displaystyle \frac{X_w-{\mu }_w}{{\sigma }_w}}\right)\hfill \end{array}$$

(15)

The sensitivities to shock parameters, $\lambda$, ${\mu }_w$, and ${\sigma }_w$, reveal how changes in these parameters impact the reliability. The magnitude and sign of these sensitivities indicate the degree and direction of the change, respectively.

2.4.2. Sensitivity of Reliability with Respect to Degradation Parameters

Similarly, for degradation parameters, such as ${\mu }_i$, ${\sigma }_i$, $b_i$, and ${\mu }_h$, ${\sigma }_h$, the sensitivities of reliability can be obtained as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\mu }_i}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right){\Lambda }_i\left(t\right)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\sigma }_i}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right){\displaystyle \frac{{\Lambda }_i\left(t\right)}{{\sigma }_i}}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial b_i}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right)\left(ln\left(t\right){\Lambda }_i\left(t\right){\left(\lambda t\right)}^i\right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\mu }_h}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right)\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\sigma }_h}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right){\displaystyle \frac{i{\Lambda }_i\left(t\right)}{\sqrt{{\sum }_{i=1}^n{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\hfill & \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial R_e\left(t\right)}{\partial {\sigma }_{\epsilon }}}& =-{\displaystyle \frac{exp(-\lambda t){\left(\lambda t\right)}^i}{i!}}\varphi \left({\displaystyle \frac{X_f-{\sum }_{i=1}^s(x_{t_i}+{\mu }_i{\Lambda }_i\left(t\right))+i{\mu }_h}{\sqrt{{\sum }_{i=1}^s{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\right){\displaystyle \frac{{\sigma }_{\epsilon }}{\sqrt{{\sum }_{i=1}^s{\sigma }_i^2{\Lambda }_i\left(t\right)+i{\sigma }_h^2+{\sigma }_{\epsilon }^2}}}\hfill \end{array}$$

(21)

Through the sensitivity analysis, we can identify which parameters play the most crucial role in system reliability. Those parameters should be the focus of the reliability improvement efforts, whether through improved design, better maintenance strategies, or improved quality control.

3. Parameter Estimation Method

The aim of this section is to achieve the parameter estimation of the proposed model. Consider that the arrival of shock affects the degradation process by not only accelerating it but also modifying the degradation pattern. Suppose that the degradation process comprises s states, with each state being initiated due to the shock arrival, lasting within the time interval $[t_s,t_{s+1}]$. If the degradation process of each state is modeled by the time transform Wiener process then the degradation path of the s-state can be expressed as:

$$X_s\left(t\right)=x_{t_s}+{\mu }_s{\Lambda }_s\left(t\right)+{\sigma }_{s}B\left({\Lambda }_s\left(t\right)\right)+ϵ,t_s\le t\le t_{s+1}$$

(22)

Here, each state s has a set of degradation observations $X_s\left(t\right)=\{x_{s,1},x_{s,2},\dots ,x_{s,n_s}\}$, with $s=0,1,2,\dots ,n$ and ${\Lambda }_s\left(t\right)=t^{b_s}$, and $s=0$ denotes the initial state, i.e., the state before the arrival of shock, and $n_s$ represents the total observations for state s prior to the arrival of the next shock. Hence, the independent increments of the Wiener process $\Delta {\mathit{X}}_s=x_{s,j+1}-x_{s,j}$, (where $j=1,\dots ,n_s-1$) follows a normal distribution with mean ${\mu }_s\Delta {\Lambda }_{s,j}$ and variance ${\sigma }_s^2\Delta {\Lambda }_{s,j}+{\sigma }_{ϵ}^2$, where $\Delta {\Lambda }_{s,j}={\Lambda }_s\left(t_{s,j+1}\right)-{\Lambda }_s\left(t_{s,j}\right)$ and the likelihood function for the observations in state s can be expressed as:

$$L_s(\Delta {\mathit{X}}_s|{\mu }_s,{\sigma }_s,{\sigma }_{ϵ})=\prod _{j=1}^{n_s-1}{\displaystyle \frac{1}{\sqrt{2\pi ({\sigma }_s^2\Delta {\Lambda }_{s,j}+{\sigma }_{ϵ}^2)}}}exp\left(-{\displaystyle \frac{{[x_{s,j+1}-x_{s,j}-{\mu }_s\Delta {\Lambda }_{s,j}]}^2}{2({\sigma }_s^2\Delta {\Lambda }_{s,j}+{\sigma }_{ϵ}^2)}}\right)$$

(23)

If the observed shock lengths $h=\{h_1,h_2.\dots ,h_s\}$ and the magnitude of shocks $W=\{W_1,W_2,\dots ,W_s\}$ are assumed to follow a normal distribution $h\sim ({\mu }_h,{\sigma }_h)$ and $W\sim ({\mu }_w,{\sigma }_w)$, respectively, then complete likelihood for the function can be expressed as:

$$\begin{array}{cc}\hfill & L(\Delta {\mathit{X}}_s,W,h,N\left(t\right)|\mathit{\mu },\mathit{\sigma },{\mu }_w,{\sigma }_w,{\mu }_h,{\sigma }_h,{\sigma }_{ϵ},\lambda ,b_i)=\hfill \\ \hfill & \prod _{i=0}^s\left[\left[\prod _{j=1}^{n_s-1}{\displaystyle \frac{1}{\sqrt{2\pi ({\sigma }_s^2\Delta {\Lambda }_{s,j}+{\sigma }_{ϵ}^2)}}}exp\left(-{\displaystyle \frac{{[x_{s,j+1}-x_{s,j}-{\mu }_s\Delta {\Lambda }_{s,j}]}^2}{2({\sigma }_s^2\Delta {\Lambda }_{s,j}+{\sigma }_{ϵ}^2)}}\right)\right]\right.\hfill \\ \hfill & \left.\times I(i>0)\left[{\displaystyle \frac{1}{2\pi {\sigma }_w{\sigma }_h}}·exp\left(-{\displaystyle \frac{{(W_i-{\mu }_w)}^2}{2{\sigma }_w^2}}-{\displaystyle \frac{{(h_i-{\mu }_h)}^2}{2{\sigma }_h^2}}\right)\right]\times \left(e^{-\lambda t}{\displaystyle \frac{{\left(\lambda t\right)}^i}{i!}}\right)\right]\hfill \end{array}$$

where $\mathit{X}=\{X_0,X_1,\dots ,X_s\}$, $\mathit{\mu }=\{{\mu }_0,{\mu }_1,\dots ,{\mu }_s\},$ $\mathit{\sigma }=\{{\sigma }_0,{\sigma }_1,\dots ,{\sigma }_s\}$, and $I(·)$ is an indicator function.

Considering the complexity of our likelihood function, which involves high- dimensional parameters, we adopted the ABC-Gibbs algorithm to estimate our model parameters. The ABC-Gibbs algorithm is a sophisticated variant of the ABC algorithm, specifically designed to handle challenges posed by complex likelihood functions with high-dimensional parameters.

The ABC-Gibbs method offers a powerful solution to estimating the parameters of complex models where the direct computation of the likelihood function is difficult or computationally intensive. This hybrid technique brings together the strengths of ABC, which circumvents the need for explicit likelihood calculation, and Gibbs sampling, a Markov-chain Monte Carlo (MCMC) method that generates samples from the joint posterior distribution of the parameters [53]. In each iteration of the ABC-Gibbs sampler, one parameter is selected for potential update. Given the current values of the other parameters, synthetic data are generated under the model. The similarity between these synthetic data and the observed data is evaluated using a pre-defined distance metric. If this distance is below a certain threshold, the proposed parameter value is accepted, otherwise it is rejected. This procedure is repeated for each parameter in turn, and over many iterations, it yields samples from the approximate joint posterior distribution of all parameters.

The ABC-Gibbs sampler allows for the estimation of the posterior distributions of each model parameter, providing detailed insights into the associated parameter uncertainties. To describe the ABC-Gibbs algorithm, let us assume that the likelihood function is obtained based on multiple states (i.e., $s=0,1,2$). The resulting ABC-Gibbs algorithm is summarized in Algorithm 1. In this algorithm, $P\left({\theta }_i\right|S\left(X_{\mathrm{obs}}\right),{\theta }_{-i})$ represents the conditional posterior distribution of the parameter ${\theta }_i$ in the s-state process, given the current values of the other parameters ${\theta }_{-i}$ and the observed data $X_{\mathrm{obs}}$. The function $S\left(X_{\mathrm{obs}}\right)$ computes the summary statistics of the observed data, which are used to compare with the synthetic data $S\left(X\right)$ generated by the model. The model $\mathcal{M}\left(\mathit{\theta }\right)$ is a function that simulates or generates synthetic data based on the current set of parameters $\mathit{\theta }$. The distance metric $\rho (S\left(X_{\mathrm{obs}}\right),S\left(X\right))$ defines the measure of dissimilarity between the summary statistics of the observed data $S\left(X_{\mathrm{obs}}\right)$ and those of the synthetic data $S\left(X\right)$. The threshold $\epsilon$ determines the maximum allowable difference between the summary statistics for accepting or rejecting the proposed parameter values. For more detail regarding ABC-Gibbs, refer to [10,50,54,55].

Algorithm 1 ABC-Gibbs algorithm for s-state degradation-shock process
Inputs:
i.	Observed dataset $X_{\mathrm{obs}}$,
ii.	Summary statistics of the observed data $S\left(X_{\mathrm{obs}}\right)$,
iii.	Desired number of samples $T>0$,
iv.	Distance measure $\rho (.)=\parallel ·\parallel$,
v.	Tolerance threshold $\epsilon \ge 0$,
vi.	Starting points ${\mathit{\theta }}^{\left(0\right)}=({\theta }_1^{\left(0\right)},\dots ,{\theta }_D^{\left(0\right)})$
Simulation:
1: for $t=1$ to T do   2:     for $i=1$ to D do   3:         Update each parameter i for the s-state conditionally on the others: $${\theta }_i^{(t+1)}\sim P\left({\theta }_i^{\left(t\right)}\right|S\left(X_{\mathrm{obs}}\right),{\theta }_1^{(t+1)},\dots ,{\theta }_{i-1}^{(t+1)},{\theta }_{i+1}^{\left(t\right)},\dots ,{\theta }_D^{\left(t\right)})$$   4:     end for   5:     Set the vector of current parameters: $${\mathit{\theta }}^{(t+1)}=({\theta }_1^{(t+1)},\dots ,{\theta }_D^{(t+1)})$$   6:     Generate synthetic data $X=\mathcal{M}\left({\mathit{\theta }}^{(t+1)}\right)$ and compute the summary statistics $S\left(X\right)$   7:     Calculate the distance $\rho (S\left(X_{\mathrm{obs}}\right),S\left(X\right))=\parallel S\left(X_{\mathrm{obs}}\right)-S\left(X\right)\parallel$   8:     if $\rho (S\left(X_{\mathrm{obs}}\right),S\left(X\right))\le \epsilon$ then   9:         Accept the parameter values 10:     else 11:         Reject the parameter values 12:     end if 13: end for

4. Illustrative Examples

In this section, we introduce both simulation and real data application examples to illustrate the effectiveness of our parameter estimation approach and the applicability of our proposed method.

4.1. Simulation Study

Consider a degradation-shock process involving four states, namely $s=0,1,2,3$. The shocks are assumed to arrive at a rate determined by a Poisson process with an intensity of $\lambda =0.004$. The degradation observations are taken every 2 h until the arrival of a shock, which marks the beginning of the next state, and the total length of time of observations is 100 h. Furthermore, both the magnitude W and the sudden jump h that arise due to a shock’s arrival are treated as independent and identically distributed (i.i.d.) normal variables. Specifically, the distribution parameters for these variables are set as ${\mu }_w=0.150,{\sigma }_w=0.020$ for W, and ${\mu }_h=0.500,{\sigma }_h=0.500$ for h. The degradation parameters of each of the four states are denoted by vectors $\mathit{\mu }=\{1.920,2.200,5.500,6.000\}$, $\mathit{\sigma }=\{0.580,0.530,2.700,2.000\}$, and $\mathit{b}=\{1.200,1.560,3.000,2.500\}$; we also incorporate noise with ${\mathit{\sigma }}_{\mathit{ϵ}}=\{0.320,0.380,0.4,0.5\}$. The process is defined as having failed if either the magnitude of a shock W exceeds $X_w=0.320$ or the degradation threshold reaches $X_f=7.400$. The objective of this simulation study is to validate our parameter estimation approach within the context of the proposed model. For this purpose, we utilize the ABC-Gibbs method, as detailed in Algorithm 1, under four sets of tolerance thresholds $\epsilon =\{0.1000,0.0500,0.0050,0.0001\}$, for $1.0\times {10}^5$ iterations. The priors of the unknown parameters are assumed to be:

• ${\mu }_w\sim N({\alpha }_{{\mu }_w},{\beta }_{{\mu }_w})$
• ${\sigma }_w\sim \mathrm{gamma}({\alpha }_{{\sigma }_w},{\beta }_{{\sigma }_w})$
• ${\mu }_h\sim N({\alpha }_{{\mu }_h},{\beta }_{{\mu }_h})$
• ${\sigma }_h\sim \mathrm{gamma}({\alpha }_{{\sigma }_h},{\beta }_{{\sigma }_h})$
• ${\mu }_s\sim N({\alpha }_{{\mu }_s},{\beta }_{{\mu }_s})$
• ${\sigma }_s\sim \mathrm{gamma}({\alpha }_{{\sigma }_s},{\beta }_{{\sigma }_s})$
• $b_s\sim N({\alpha }_{b_s},{\beta }_{b_s})$
• ${\sigma }_{ϵ}\sim \mathrm{gamma}({\alpha }_{{ϵ}_s},{\beta }_{{ϵ}_s})$

where ${\alpha }_{{\mu }_w},{\beta }_{{\mu }_w},{\alpha }_{{\sigma }_w},{\beta }_{{\sigma }_w},{\alpha }_{{\mu }_h},{\beta }_{{\mu }_h},{\alpha }_{{\sigma }_h},{\beta }_{{\sigma }_h},{\alpha }_{{\mu }_s},{\beta }_{{\mu }_s},{\alpha }_{{\sigma }_s},{\beta }_{{\sigma }_s},{\alpha }_{b_s},{\beta }_{b_s}$, and ${\alpha }_{{ϵ}_s},{\beta }_{{ϵ}_s}$ are known hyper-parameters. The simulation results for both shock and degradation parameters are presented in

Table 1. Summary of the estimated mean of each shock parameter sample generated using ABC-Gibbs of the simulated datasets.

			Sample Mean
	No.Itr	$\mathit{\epsilon }$	$\mathit{\lambda }$	${\mathit{\mu }}_{\mathit{h}}$	${\mathit{\sigma }}_{\mathit{h}}$	${\mathit{\mu }}_{\mathit{w}}$	${\mathit{\sigma }}_{\mathit{w}}$
TVs	$1.0\times {10}^5$	-	4.0000	0.5000	2.0000	0.1500	0.0200
		0.1000	3.9720	0.5311	2.0426	0.1679	0.0226
		0.0500	3.9717	0.5317	2.0426	0.1652	0.0223
		0.0050	3.9760	0.5307	2.0430	0.1642	0.0223
		0.0010	3.9625	0.5123	2.0428	0.1595	0.0223

TVs = true values. No.Itrs = number of iterations.

and

Table 2. Summary of the estimated mean of each degradation parameter sample generated using ABC-Gibbs of the simulated datasets.

			Sample Mean
			State0				State1
	No.Itrs	$\mathit{\epsilon }$	${\mathit{\mu }}_{\mathbf{0}}$	${\mathit{\sigma }}_{\mathbf{0}}$	b0	${\mathit{\sigma }}_{\mathit{ϵ}\mathbf{0}}$	${\mathit{\mu }}_{\mathbf{1}}$	${\mathit{\sigma }}_{\mathbf{1}}$	b1	${\mathit{\sigma }}_{\mathit{ϵ}\mathbf{1}}$
TVs	$1.0\times {10}^5$	-	1.9200	0.5800	1.2000	0.3200	2.2000	0.5300	1.5600	0.3800
		0.1000	1.9724	0.5857	1.2550	0.3056	2.1658	0.5211	1.5931	0.3857
		0.0500	1.9754	0.5853	1.2513	0.3052	2.1685	0.5217	1.5912	0.3855
		0.0050	1.9558	0.5778	1.2428	0.3074	2.1649	0.5231	1.5968	0.3858
		0.0001	1.9128	0.5815	1.2154	0.3125	2.1892	0.5282	1.5681	0.3836
			Sample Mean
			State2				State3
	No.Itrs	$\mathit{\epsilon }$	${\mathit{\mu }}_{\mathbf{2}}$	${\mathit{\sigma }}_{\mathbf{2}}$	b2	${\mathit{\sigma }}_{{\mathit{ϵ}}_{\mathbf{2}}}$	${\mathit{\mu }}_{\mathbf{3}}$	${\mathit{\sigma }}_{\mathbf{3}}$	b1	${\mathit{\sigma }}_{\mathit{ϵ}\mathbf{3}}$
TVs	$1.0\times {10}^5$	-	5.5000	2.7000	0.3000	0.4000	6.0000	2.0000	2.5000	0.5000
		0.1000	5.6195	2.7521	0.3799	0.3916	5.8990	1.9596	2.8011	0.5120
		0.0500	5.6239	2.7510	0.3755	0.3915	5.8988	1.9626	2.8065	0.5119
		0.0050	5.6298	2.7580	0.3525	0.3914	5.8630	1.9892	2.6846	0.5133
		0.0001	5.5679	2.7421	0.3104	0.3945	5.9845	2.0126	2.5110	0.5082

. The plots of the generated samples drawn for each parameter with respect to the summary statistics and the ABC-Gibbs posterior density of each parameter are depicted in Figure 3 (for shock) and Figure 4, Figure 5, Figure 6 and Figure 7 (for degradation).

Upon analyzing the results in Table 1 and Table 2, it becomes evident that the parameter estimates derived through the ABC-Gibbs method closely approximate the true values. Furthermore, the accuracy of these estimates is significantly impacted by the tolerance threshold selected within the ABC-Gibbs algorithm. Notably, as the tolerance threshold is lowered, the estimated values increasingly align with the true values. This highlights the critical importance of choosing an appropriate tolerance threshold for achieving precise parameter estimates.

The scatter plots depicted in Figure 3, Figure 4, Figure 5 and Figure 6 offer a visual representation of the relationship between the summary statistics and the samples drawn from the conditional posterior distribution of each parameter. These plots distinguish acceptance and rejection regions with different colors, as described in the legend. The acceptance region comprises summary statistics that fulfill the tolerance criteria set by the ABC-Gibbs algorithm, while the rejection region includes those that do not meet these criteria. A visual inspection of these figures reveals a clear separation between acceptance and rejection regions, indicating that the selected summary statistics effectively capture the underlying data-generation process. The clustering of samples around the accepted summary statistics further supports the validity and effectiveness of these statistics in representing the key characteristics of the data. Similarly, the density plots in Figure 3, Figure 4, Figure 5 and Figure 6 provide a succinct visualization of parameter distributions using samples generated by the ABC-Gibbs algorithm across various tolerance thresholds. These plots illustrate a trend where the densities of each parameter converge toward a normal distribution, consistent with the central limit theorem. According to this theorem, the sum of a large number of independent and identically distributed (i.i.d.) random variables tends to converge to a normal distribution, regardless of the initial distribution of the individual variables [56].

4.2. Practical Application

A practical implementation of MEMS, which was performed at Sandia National Laboratories, is commonly referenced to demonstrate degradation-shock models in the literature [37,57,58]. Given its prevalent use and illustrative value, this section will also employ MEMS as a case study to showcase and discuss the effectiveness and applicability of the proposed model.

MEMS micro-engines are primarily susceptible to two competing failure mechanisms: (1) soft failure, primarily due to continuous wear and debris accumulation, and (2) hard failure, typically caused by hub fractures resulting from random shocks. As the micro-engine degrades, its resistance to random shocks decreases, and wear accelerates, potentially altering the degradation pattern. In light of this understanding, we apply our proposed reliability model, which accounts for changes in the degradation pattern. This model offers an alternative approach for predicting the micro-engine’s lifespan and operational reliability. All parameter values and their sources are listed in

Table 3. All parameter values and their sources.

Parameters	Values	Sources
$X_f$	0.00125 μm3	Tanner et al. [57]
$X_w$	1.5 Gpa	Rafiee et al. [59]
${\mu }_0$	$8.4823\times {10}^{-9}$ μm3	Peng et al. [60]
${\sigma }_0$	$6.00\times {10}^{-10}$ μm3	Peng et al. [60]
${\mu }_1$	$10.4823\times {10}^{-9}$ μm3	Assumption
${\sigma }_1$	$3.00\times {10}^{-10}$ μm3	Assumption
${\mu }_2$	$12.4823\times {10}^{-9}$ μm3	Assumption
${\sigma }_2$	$4.00\times {10}^{-10}$ μm3	Assumption
$x_{t1}$	0.00021 μm3	Assumption
$x_{t2}$	0.00063 μm3	Assumption
${\mu }_h$	$1.2\times {10}^{-4}$ μm3	An and Sun [61]
${\sigma }_h$	$4.0\times {10}^{-5}$ μm3	An and Sun [61]
${\mu }_w$	1.2 Gpa	An and Sun [61]
${\sigma }_w$	0.4 Gpa	An and Sun [61]
${\sigma }_{ϵ}$	0.14 Gpa	Assumption
$\lambda$	$5\times {10}^{-5}$	Tanner and Dugger [62]
$b_1=b_2=b_3$	1	Assumption

. It is important to note that although the parameter values are derived from multiple complex studies, the failure dataset used is from a single MEMS experiment. The analysis results are presented in Figure 8, Figure 9, Figure 10 and Figure 11.

It can be observed that the reliability function $Re\left(t\right)$ depicted in Figure 8 demonstrates the system’s decreasing reliability over time. Initially, $Re\left(t\right)$ is close to 1, indicating a high likelihood of the system being functional. However, as time progresses, the reliability function decreases monotonically, eventually approaching zero. This behavior is expected due to the cumulative effects of both the degradation process and random shocks that progressively deteriorate the system’s performance. The figure shows a notable change in the slope of the reliability function; this inflection point suggests a change in the degradation pattern, possibly indicating a shift from a predominantly slow degradation phase to a phase where degradation accelerates or the influence of random shocks becomes more pronounced. This change in pattern highlights the importance of considering both gradual degradation processes and sudden shocks in reliability modeling. Figure 9 and Figure 10 illustrate the sensitivity of the reliability function $Re\left(t\right)$ with respect to the degradation-shock parameters across different states. The results reveal that the sensitivity of reliability increases progressively with each state transition. Specifically, ${\mu }_0$ and ${\sigma }_0$ exhibit lower sensitivity compared to ${\mu }_1$ and ${\sigma }_1$, which in turn are less sensitive than ${\mu }_2$ and ${\sigma }_2$. This trend underscores the impact of shocks on the degradation process: the first state, prior to any shocks, follows a relatively linear and predictable degradation path, resulting in lower sensitivity to parameter changes. Following the first shock, the system enters the second state, where the degradation path becomes more variable and non-linear, increasing sensitivity to changes in ${\mu }_1$ and ${\sigma }_1$. The second shock further exacerbates the degradation process, leading to the third state, where the degradation becomes highly non-linear and complex, making the system most sensitive to ${\mu }_2$ and ${\sigma }_2$. Additionally, the parameters $b_1$, $b_2$, and $b_3$ represent the time transformation factors in the degradation process, indicating the degree of non-linearity in each state. The analysis shows that $b_3$ is significantly more sensitive than $b_1$ and $b_2$, highlighting the increased non-linear effects in the third state. This increased sensitivity of $b_3$ reflects the compounded impact of two shocks, which introduce substantial deviations from the initial degradation pattern and result in a more complex degradation trajectory.

The shock parameters $\lambda$, ${\mu }_w$, and ${\sigma }_w$ also play a crucial role in influencing the reliability of the system. The sensitivity analysis indicates that these shock parameters significantly affect the reliability, particularly in the later states. As the rate of shock occurrence ($\lambda$) and the magnitude of the shocks (${\mu }_w$ and ${\sigma }_w$) increase, the degradation process becomes more abrupt and unpredictable, thereby increasing the sensitivity of the reliability function. The impact of these shock parameters is more pronounced in the second and third states, where the system has already undergone one or more shocks, leading to a heightened response to further shocks.

We further investigate the effect of $b_1$, $b_2$, and $b_3$ within the reliability curve to understand how these time transformation factors influence system reliability across different states. The result is presented in Figure 11. The analysis reveals that $b_1$ does not significantly affect the reliability function in the first state, indicating that initial degradation processes are relatively insensitive to variations in $b_1$ within the selected values. In contrast, the factors $b_2$ and $b_3$ exhibit a pronounced impact on the reliability function in the second and third states, respectively. As $b_2$ increases, the reliability function declines more rapidly, suggesting that post-shock degradation in the second state is highly sensitive to changes in $b_2$. Similarly, increasing $b_3$ results in a faster decrease in reliability, reflecting the significant acceleration in degradation after the second shock. This accelerated degradation after the second shock can be attributed to the accumulation of shock impacts, which compound the damage and introduce more severe degradation patterns.

These observations align with the theoretical expectations that shocks not only accelerate degradation but also introduce significant changes in the degradation patterns, thereby increasing the system’s sensitivity to parameter variations as it transitions through successive states. This comprehensive sensitivity analysis, encompassing both degradation and shock parameters, is crucial for understanding the robustness of the reliability model and for identifying key parameters that influence system performance, thereby guiding targeted reliability improvement efforts.

5. Conclusions

This study presents a comprehensive multi-state reliability model for systems subject to competing failure processes driven by both degradation and random shocks. By considering state-varying degradation patterns, our model captures the dynamic and complex nature of real-world degradation processes. The effectiveness of our parameter estimation approach was validated through a simulation study. This simulation involved a complex system influenced by a degradation-shock process with multiple states. The shocks, modeled by a Poisson process, and the degradation parameters were tested under various tolerance thresholds using the ABC-Gibbs method. The simulation results demonstrated that the parameter estimates closely approximate the true values, emphasizing the importance of selecting appropriate tolerance thresholds for achieving accurate results. In addition to the simulation, a practical application was performed using a MEMS case study, which illustrated the applicability and effectiveness of the proposed model. The reliability function analysis revealed that the system’s reliability decreases over time, with notable changes in the slope indicating shifts in the degradation pattern. Sensitivity analysis further highlighted the increasing sensitivity of the reliability function to parameter changes across different states, underscoring the compounded impact of shocks and the non-linear degradation trajectory in the later states.

Overall, the proposed multi-state reliability model, combined with sensitivity analysis and parameter estimation using the ABC-Gibbs method, offers a robust framework for analyzing complex systems with degradation-shock dependencies. This approach not only enhances the accuracy of reliability predictions but also guides targeted reliability improvement efforts, making it a valuable tool for engineers and researchers in the field of reliability engineering. Furthermore, the model’s adaptability has significant potential across various disciplines. In materials science, it can be used to predict material degradation and failure, while in operational research, it could optimize resource allocation and maintenance strategies for systems subject to varying degradation patterns. Future work will focus on extending the model to incorporate additional failure mechanisms, exploring its applicability to a broader range of engineering systems, and refining the parameter estimation techniques to handle even more complex and high-dimensional scenarios. Moreover, integrating machine learning algorithms with our model could further enhance predictive capabilities and facilitate real-time reliability assessments, ultimately contributing to more resilient and reliable engineering products across diverse industries such as aerospace, automotive, and heavy machinery.