An artificial neural network supported Wiener process based reliability estimation method considering individual difference and measurement error

Abstract

Due to the powerful ability of artificial neural network in data fitting, it has been applied to describe the mean function in Wiener process for degradation modeling and estimating reliability. However, the previously published method neglects the measurement error and individual difference, which need to be considered in reliability estimation. In order to handle the above issues, the artificial neural network supported Wiener process is improved. Both the individual difference and measurement error are considered. The measurement error is described by a normal distribution with zero mean. The individual difference is described by introducing hyper parameters. In order to demonstrate the proposed method, a simulation study and a case study have been performed. It can be seen that, considering the measurement error and individual difference can improve the reliability and lifetime estimation accuracies, resulting in the proposed method is more powerful and more practical in engineering practice.

Introduction

Nowadays, with the development of manufacturing technique, modern products always have long lifetime and high reliability [1,2]. For some safety-critical components, the high reliability during the whole lifetime is demanded to guarantee the safety of the system [3]. Hence, before the components are put into the system the reliability tests are usually performed [4]. Using the reliability tests, two kinds of data can be obtained, including degradation data and lifetime data [5]. In engineering practice, due to the long lifetime of modern products, degradation data is more easily to be collected than lifetime data by reliability testing. Hence, degradation modeling based methods attract lots of attentions in reliability estimation [6].

In degradation modeling based reliability estimation, stochastic process is one of the most practical approaches [7]. There are three kinds of widely used stochastic processes. Wiener process can be used to model both monotonic and non-monotonic degradation processes, while Gamma process and inverse Gaussian process are only suitable for modeling monotonic degradation processes [8], [9], [10], [11]. Three factors need to be determined when applying stochastic process, including degradation process, mean function and model parameters [12].

In most cases, there may be several suitable processes for a specific degradation dataset [13]. Several methods have been presented to select and combine the candidate processes. Several criterions have been proposed to assess the candidate processes and select the best process, such as Akaikes information criterion (AIC), Bayesian information criterion (BIC), etc. [14], [15], [16]. Furthermore, several model combining methods have also been presented to combine the candidate processes, such as adjustment factor method, Bayesian model averaging method (BMA) and evidence theory based method [17], [18], [19].

Bayesian method has been applied to estimate the model parameters, due to its powerful ability in data fusion and parameters estimation [20], [21], [22], [23]. For example, Ref. [23] has presented several inverse Gaussian process models and the related Bayesian analysis methods to model the degradation process and infer the degradation of the product. Besides the Bayesian approaches, some other methods have also been proposed to adapt stochastic process to small sample conditions, such as random fuzzy theory based methods, evidence theory based methods, etc. [24,25].

For a specific degradation dataset, there may be a number of available mean functions, such as multi-phase based mean functions, mean functions with time varying degradation rates, etc. [26], [27], [28], [29]. When using stochastic process, the mean function is normally determined by model selecting method or combination method among the existing candidates. However, the existing mean functions should not be suitable for all applications. For a specific application, it may be not suitable to determine the mean function by selecting or combining within the existing candidates. Hence, the mean functions are determined by degradation mechanism analysis for some specific applications. Ma et al. have presented an engineering‐driven performance degradation analysis method [30]. The degradation of a kind of piston pumps is modeled by inverse Gaussian process. The mean function is given by analyzing the degradation mechanism of the pump. Similar methods can also be found in Ref. [31]. The researchers argue that the theoretical mean functions predicted by analyzing the degradation mechanism are not the existing candidates and the theoretical mean functions show superiority in engineering practice. However, it is obvious that the mechanism analysis based methods need large amount of knowledge about the application objects.

Due to the flexibility on relationship describing and fitting, artificial neural network (ANN) has been widely used and shows superiority in some practical applications. For instance, Hong et al. have employed ANN in oil debris detection and proposed a hybrid detection strategy to reduce the detecting error [32]. An experimental study is carried out to verify the proposed method. It is concluded that the detection accuracy is improved due to utilizing the advantage of ANN. Furthermore, ANN has also been introduced to describe the relationship between the wear rate and its affecting factors in wear failure analysis [33,34]. It is concluded that the ANN supported models can predict the wear failure process precisely. Ma et al. have constructed a degradation modeling approach for proton exchange membrane fuel cell based on deep learning method [35]. The recurrent neural network with grid long short-term memory cell is utilized in the degradation prognostic model. The superiorities of the proposed method are verified based on an actual degradation dataset from eight different fuel cells with eight different load profiles. The above mentioned researches have validated the strong ability of ANN on relationship description and data fitting. Hence, in order to describe the mean function, we have combined ANN with Wiener process for estimating reliability of the products [36]. ANN is applied to fit the mean function in Wiener process. Simulation studies have been performed to verify the advantages of the presented method. In addition, an actual case is also used to demonstrate the presented reliability estimation method. It is concluded that ANN is capable in mean function fitting and the presented method is capable in reliability estimation. However, the measurement error and individual difference are neglected in the above mentioned research.

In fact, the measurement error is an inevitable phenomenon and need to be considered. Ma et al. have proposed an inverse Gaussian process model with measurement error and the corresponding parameters estimation method based on expectation maximization algorithm [37]. The presented model and the associated method are verified by comparing to the method neglecting measurement error based on a simulation study and an actual case study. The authors argue that considering the measurement error is necessary in stochastic process based reliability estimation. Ye et al. have proposed several Wiener process models with measurement errors and the corresponding statistical inference methods to capture the effects of the measurement errors on the degradation modeling and reliability estimation [38]. The necessity of considering the measurement errors is verified based on a simulation study and two illustrative examples. In Ref. [39], the imperfect inspection issue in dynamic reliability assessment is focused and handled, especially the measurement error. The dynamic reliability assessment for nonrepairable multistate systems is discussed. A set of two-stage recursive Bayesian formulations are constructed to utilize the imperfect inspection data. The effectiveness of the proposed method is verify based on a case study. Li et al. have presented a Wiener process model considering the measurement errors to analyze accelerated degradation data [40]. The measurement error term is introduced in Wiener process model and described by a normal distribution with zero mean. The maximum likelihood estimations of the model parameters are deduced. The superiority of the proposed method is illustrated by comparing to the method neglecting the measurement errors based on an actual degradation dataset regarding LED. Besides, Pan et al. have adopted expectation maximization algorithm to the Wiener process degradation model with measurement errors [41]. It is argued that, the effectiveness of the Wiener process based method has been improved due to considering the measurement error. In Ref. [42], inverse Gaussian process model has been extended to consider the measurement errors by introducing measurement error item, which is assumed to follow a normal distribution with zero mean. The model parameters are estimated based on Monte Carlo algorithm. A simulation study and two actual case studies are performed to verify the effectiveness of the proposed method. In brief, the measurement error is an inevitable issue and need to be considered in degradation modeling and reliability estimation.

Besides the measurement error, the individual difference also need to be considered. In engineering practice, hyper parameters are always introduced to handle the random effects caused by individual difference. In Ref. [43], several Gaussian processes with random effects have been introduced. The unit-to-unit variability is incorporated into the model by applying the natural conjugate distribution. The presented models are demonstrated based on a Monte Carlo simulation study and several case studies. Based on the simulation results, the author argues that the presented models provide the applicability and effectiveness in degradation modeling, due to considering the individual difference. Wang has presented a class of Wiener processes with random effects to analyze degradation data [44]. The model parameters can vary from unit to unit by introducing the random drift and diffusion parameters. Expectation maximization algorithm is applied to estimate the unknown parameters. Based on the simulation study, it is concluded that introducing random effects in Wiener processes is necessary. In addition, a series of gamma process models with the random effects have also been presented in Ref. [45]. The drift and diffusion parameters are assumed to be randomly distributed to describe the random effects caused by the individual difference. The model parameters are estimated by Bayesian inference. The necessity of considering the individual differences is proved by several case studies. Based on the above discussions, it can be seen that considering the individual difference is necessary.

As discussed above, in degradation modeling and reliability estimation, both the measurement error and individual difference should to be considered. Hence, the above two issues are focused in the presented research. Furthermore, because of the modeling ability on both non-monotonic and monotonic degradation processes, Wiener process is presented as an example in this paper. Hence, the ANN supported Wiener process based reliability estimation method is improved. The distinguished features of the proposed reliability estimation method are as follows. Both the measurement error and individual difference are considered in ANN supported Wiener process based degradation model. The measurement error is described by a normal distribution with zero mean. The individual difference is described by introducing hyper parameters. The ANN supported mean function is trained by maximizing the likelihood of the corresponding Wiener process model to the entire evaluation degradation dataset. The hyper parameters for the testing sample population and the model parameters for the monitoring individual are evaluated by Bayesian inference combining with maximum likelihood estimation (MLE).

This paper is organized as follows. The ANN supported Wiener process considering the measurement error and individual difference is presented in Section 2. A three stages parameter estimation method and the associated procedure are presented in Section 3. The proposed model and the corresponding reliability estimation method are demonstrated based on a simulation study and an illustrate example in Sections 4 and 5. Finally, the conclusions are presented in Section 6.

For a Wiener degradation process Y(t), Y(t) has independent increments. Furthermore, the increments ΔY follow Gaussian distribution, as

$$\Delta Y\left(t\right)\sim N\left(\Delta \Lambda \left(t\right),{\lambda }^2\Delta \Lambda \left(t\right)\right)$$

where N(•,•) indicates Gaussian distribution. ΔΛ(t)=Λ(t+Δt)-Λ(t), where Δt is time increment. Λ(t) is the mean function and needs to be monotonically increasing.

Furthermore, the degradation at time t is

$$Y\left(t\right)\sim N\left(\Lambda \left(t\right),{\lambda }^2\Lambda \left(t\right)\right)$$

where Λ(t) and λ²Λ(t) means the degradation mean and variance at time t. λ is called diffusion parameter. The parameter vector for Wiener process θ_W = (λ).

An ANN supported Wiener process with measurement error and individual difference is proposed. Considering the measurement error, the degradation observation z is given by Eq. (3).

$$z=y+\epsilon$$

where y means the true degradation. Measurement error term ε is assumed to follow a Gaussian distribution with σ² variance and zero mean. The measurement error parameter vector θ_E=(σ). Furthermore, similarly as the assumptions in Refs. [37,38,[40], [41], [42]], the measurement error only depend on the utilized sensor and irrelevant to degradation process, so the measurement error is assumed to be s-independent with the other process model's parameters.

Based on an ANN and two parameters, the mean function is given as

$$\Lambda \left(t\right)=\mu ANN\left(\alpha t\right)$$

where ANN(•) means a single input single output ANN, which is used to indicate the failure mechanism and does not vary with different individuals. μ and α are introduced to control the changing rate and the shape of the mean function. Mean function parameter vector θ_Λ=(α,μ), which are varied with different individuals. The schematic of the proposed ANN supported mean function is shown in Fig. 1.

In the above ANN, sigmoid function is used as activation function for the hidden neurons. Namely, the ANN can be described as

$$\{\begin{array}{c}{x}_{\mathrm{o}}=\sum _{k=1}^n{w}_k{h}_k\hfill \\ h_k=S\left(v_k{x}_{\mathrm{i}}+b_k\right)\hfill \\ S\left(x\right)=\frac{1}{1+e^{-x}}\hfill \end{array}$$

where x_o is the output of the ANN, x_i is the input of the ANN, h_k is the value of kth hidden neuron, v_k is kth weight of the input layer, w_k is kth weight of the output layer, b_k is the bias of kth hidden neuron and n is the neuron size of the hidden layer. The ANN parameter vector θ_ANN=(V,W,B) and does not vary with different individuals, where V=(v₁,v₂,…,v_n), W=(w₁,w₂,…,w_n) and B=(b₁,b₂,…,b_n). Here, the model parameter vector θ=(θ_Λ,θ_W,θ_E)=(α,μ,λ,σ). In order to consider the individual differences among the units, the model parameters are described by truncated normal (TN) distribution with hyper parameters, including the diffusion parameter λ, measurement error parameter σ, mean function parameters α and μ. Namely,

$$\{\begin{array}{c}\alpha \sim TN\left({\gamma }_{\alpha },{\delta }_{\alpha }^2,0,+\infty \right),\mu \sim TN\left({\gamma }_{\mu },{\delta }_{\mu }^2,0,+\infty \right)\hfill \\ \lambda \sim TN\left({\gamma }_{\lambda },{\delta }_{\lambda }^2,0,+\infty \right),\sigma \sim TN\left({\gamma }_{\sigma },{\delta }_{\sigma }^2,0,+\infty \right)\hfill \end{array}$$

where θ^H=(γ_α,δ_α,γ_μ,δ_μ,γ_λ,δ_λ,γ_σ,δ_σ) means the corresponding hyper parameter vector. TN(•,•,0,+∞) means the TN distribution and the corresponding probability density function (PDF) is given by

$$f_{\text{TN}}\left(x|{\gamma }_x,{\delta }_x^2,0,+\infty \right)=\{\begin{array}{c}\frac{\frac{1}{{\delta }_x}\varphi \left(\frac{x-{\gamma }_x}{{\delta }_x}\right)}{\Phi \left(\frac{{\gamma }_x}{{\delta }_x}\right)}0\le x\hfill \\ 0\text{other}\hfill \end{array}$$

where ϕ(•) and Φ(•) mean the PDF and the cumulative probability function (CDF) of standard normal distribution.

In this section, a three-stages parameter estimation method is constructed for applying the proposed ANN supported Wiener process. First, an ANN supported mean function training approach is constructed based on the likelihood function, see Stage I. Second, the priors of the hyper parameters are set based on the training results and MLE, see Stage II. Finally, Bayesian inference method is introduced to estimate the hyper parameters when evaluating the testing sample population and update the model parameters when predicting the individual, see Stage II.

For ith sample, the degradation dataset includes the time and observation vectors, Eqs. (8) and (9).

$$T_i=\left[t_{i1},t_{i2},\dots t_{i{n}_i}\right]$$

$$Z_i=\left[z_{i1},z_{i2},\dots z_{i{n}_i}\right]$$

where n_i is the measurement size of ith sample, t_ij means j-th inspection time of ith sample and z_ij is the corresponding measured degradation.

Measurement error is considered in the presented research, so the measured degradation is sum of the true degradation and the measurement error. Namely, z_ij=y_ij+ε_ij. Hence, for ith sample, the corresponding PDF is given by

$$f\left(Z_i|{\theta }_i,ANN,T_i\right)={\int }_{-\infty }^{+\infty }\cdots {\int }_{-\infty }^{+\infty }\prod _{j=1}^{n_i}{f}_{\mathrm{N}}\left(\Delta z_{ij}-\Delta {\epsilon }_{ij}|\Delta {\Lambda }_{ij},{\lambda }_i^2\Delta {\Lambda }_{ij}\right)\phi \left(\Delta {\epsilon }_i\right)d\Delta {\epsilon }_{i1}\cdots d\Delta {\epsilon }_{i{n}_i}$$

where f_N(•|•,•) means the PDF of Gaussian distribution. ΔΛ_ij=Λ(t_ij)-Λ(t_i(j-1)), Δz_ij=z_ij-z_i(j-1), Δε_ij=ε_ij-ε_i(j-1), z_i0=0 and t_i0=0. θ_i=(α_i,μ_i,λ_i,σ_i) means the model parameters for ith sample. Furthermore, Δε_i follows a multivariate normal distribution, as

$$\phi \left(\Delta {\epsilon }_i\right)=\frac{1}{{\left(2\pi \right)}^{n_i/2}{\left|\Sigma \right|}^{1/2}}\exp \left[-\frac{1}{2}{\left(\Delta {\epsilon }_i\right)}^{\mathrm{T}}{\Sigma }^{-1}\left(\Delta {\epsilon }_i\right)\right]$$

where the tridiagonal matrix Σ is given by Eq. (12) and |Σ| means the corresponding determinant covariance.

$${\Sigma }_{hl}=\{\begin{array}{c}{\sigma }^{2}h=l=1\hfill \\ 2{\sigma }^{2}h=l>1\hfill \\ -{{\sigma }^2}^{\left|h-l\right|}=1\hfill \\ 0\text{other}\hfill \end{array}$$

For ith sample, the corresponding likelihood function is given by

$$L_{\mathrm{I}}\left(Z_i|{\theta }_i,ANN,T_i\right)=f\left(Z_i|{\theta }_i,ANN,T_i\right)$$

The likelihood function for the entire degradation testing dataset is given by

$$L_{\mathrm{P}}\left(Z_{1:n_{\mathrm{d}}}|{\theta }_{1:n_{\mathrm{d}}},ANN,T_{1:n_{\mathrm{d}}}\right)=\prod _{i=1}^{n_{\mathrm{d}}}{L}_{\mathrm{I}}\left(Z_i|{\theta }_i,ANN,T_i\right)=\prod _{i=1}^{n_{\mathrm{d}}}f\left(Z_i|{\theta }_i,ANN,T_i\right)$$

where $T_{1:n_{\mathrm{d}}}=\left(T_1;T_2;...;T_{n_{\mathrm{d}}}\right)$ means the inspection times of all the n_d samples and $Z_{1:n_{\mathrm{d}}}=\left(Z_1;Z_2;...;Z_{n_{\mathrm{d}}}\right)$ means the corresponding degradation observations. ${\theta }_{1:n_{\mathrm{d}}}=\left({\theta }_1;{\theta }_2;...;{\theta }_{n_{\mathrm{d}}}\right)$ is the model parameters for all the n_d samples.

Hence, the corresponding minus log-likelihood function is given by

$$l_{\mathrm{P}}\left(Z_{1:n_{\mathrm{d}}}|{\theta }_{1:n_{\mathrm{d}}},ANN,T_{1:n_{\mathrm{d}}}\right)=-\ln L_{\mathrm{P}}\left(Z_{1:n_{\mathrm{d}}}|{\theta }_{1:n_{\mathrm{d}}},ANN,T_{1:n_{\mathrm{d}}}\right)=-\sum _{i=1}^{n_{\mathrm{d}}}\ln f\left(Z_i|{\theta }_i,ANN,T_i\right)$$

It is widely known that, the goodness-of-fit of the model to the dataset is better when the corresponding likelihood is bigger. In order words, the corresponding minus log-likelihood is smaller the goodness-of-fit is better. Hence, the loss function is constructed based on the minus log-likelihood function and a regularization term is introduced to avoid the overfitting problem. The loss function is given by

$$E=l_{\mathrm{P}}\left(Z_{1:n_{\mathrm{d}}}|{\theta }_{1:n_{\mathrm{d}}},ANN,T_{1:n_{\mathrm{d}}}\right)+K\left(\sum _{i=1}^n\left|b_i\right|+\sum _{i=1}^n\left|v_i\right|+\sum _{i=1}^n\left|w_i\right|\right)$$

where K means the regularization coefficient.

Furthermore, the mean function should be monotonically increasing. The parameters μ, α, λ and σ should be positive considering their meanings. Hence, the ANN and ${\theta }_{1:n_{\mathrm{d}}}$ can be obtained by minimizing the loss function with the above inequality constraints based on the following method, as

$$\begin{array}{c}\min E=l_{\mathrm{P}}\left(Z_{1:n_{\mathrm{d}}}|{\theta }_{1:n_{\mathrm{d}}},ANN,T_{1:n_{\mathrm{d}}}\right)+K\left(\sum _{i=1}^n\left|b_i\right|+\sum _{i=1}^n\left|v_i\right|+\sum _{i=1}^n\left|w_i\right|\right)\hfill \\ \mathrm{s}.\mathrm{t}.\frac{dANN\left(x_{\mathrm{i}}\right)}{d{x}_{\mathrm{i}}}>0\hfill \\ \forall i,{\mu }_i>0,{\alpha }_i>0,{\lambda }_i>0,{\sigma }_i>0\hfill \end{array}$$

where

$$\frac{\partial ANN\left(x_{\mathrm{i}}\right)}{\partial x_{\mathrm{i}}}=\sum _{k=1}^n{v}_k{w}_k{S}^{\prime }\left(v_k{x}_i+b_k\right)$$

where the derivative of the activation function is given by

$$S^{\prime }\left(x\right)=\frac{e^{-x}}{{\left(1+e^{-x}\right)}^2}$$

In the presented research, the above optimization problem is handled based on genetic algorithm by using ga function in MATLAB. For the simulation and case studies presented in this paper, the main configuration parameters are set as follows. The crossover fraction is 0.8. The constraint tolerance is 0.001. The population size is 200. The function tolerance is 0.000001. The trained results include the trained ANN parameter vector ${\widehat{\theta }}_{\text{ANN}}=\left(\widehat{V},\widehat{W},\widehat{B}\right)$ and trained model parameter vector ${\widehat{\theta }}_{1:n_{\mathrm{d}}}=\left({\widehat{\theta }}_1;{\widehat{\theta }}_2;...;{\widehat{\theta }}_{n_{\mathrm{d}}}\right)$, where ${\widehat{\theta }}_i=\left({\widehat{\alpha }}_i,{\widehat{\mu }}_i,{\widehat{\lambda }}_i,{\widehat{\sigma }}_i\right)$. In the following sections, $\widetilde{ANN}$ is used to indicate the trained ANN. Furthermore, the ANN only depends on the failure-generation mechanism, so it does not vary with different individuals and determined at Stage I.

Hyper parameters are introduced to describe the model parameters to consider the individual differences among the units. In the presented research, Bayesian inference is applied to estimate the hyper parameters based on the entire degradation testing dataset and the training results of the model parameters, which are used to set the prior distributions of the hyper parameters. In this section, the priors of hyper parameters are set by MLE based on the above training results of the model parameters. The MLEs of the hyper parameters are obtained by solving Eqs. (20)–(23). The corresponding deductions can be found in the supplementary material.

$$\{\begin{array}{c}\frac{1}{\Phi \left(\frac{{\gamma }_{\alpha }}{{\delta }_{\alpha }}\right)}\varphi \left(\frac{{\gamma }_{\alpha }}{{\delta }_{\alpha }}\right)\frac{{\gamma }_{\alpha }}{\sqrt{2\pi }}\exp \left(-\frac{{\gamma }_{\alpha }^2}{2{\delta }_{\alpha }^2}\right)+{\gamma }_{\alpha }=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\alpha }}_i\hfill \\ 1-\frac{1}{\Phi \left(\frac{{\gamma }_{\alpha }}{{\delta }_{\alpha }}\right)}\varphi \left(\frac{{\gamma }_{\alpha }}{{\delta }_{\alpha }}\right)\frac{{\gamma }_{\alpha }^2}{\sqrt{2\pi }{\delta }_{\alpha }^2}\exp \left(-\frac{{\gamma }_{\alpha }^2}{2{\delta }_{\alpha }^2}\right)-{\gamma }_{\alpha }^2=\frac{1}{n_d}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\alpha }}_i^2-\frac{2{\gamma }_{\alpha }}{n_d}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\alpha }}_i\hfill \end{array}$$

$$\{\begin{array}{c}\frac{1}{\Phi \left(\frac{{\gamma }_{\mu }}{{\delta }_{\mu }}\right)}\varphi \left(\frac{{\gamma }_{\mu }}{{\delta }_{\mu }}\right)\frac{{\gamma }_{\mu }}{\sqrt{2\pi }}\exp \left(-\frac{{\gamma }_{\mu }^2}{2{\delta }_{\mu }^2}\right)+{\gamma }_{\mu }=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\mu }}_i\hfill \\ 1-\frac{1}{\Phi \left(\frac{{\gamma }_{\mu }}{{\delta }_{\mu }}\right)}\varphi \left(\frac{{\gamma }_{\mu }}{{\delta }_{\mu }}\right)\frac{{\gamma }_{\mu }^2}{\sqrt{2\pi }{\delta }_{\mu }^2}\exp \left(-\frac{{\gamma }_{\mu }^2}{2{\delta }_{\mu }^2}\right)-{\gamma }_{\mu }^2=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\mu }}_i^2-\frac{2{\gamma }_{\mu }}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\mu }}_i\hfill \end{array}$$

$$\{\begin{array}{c}\frac{1}{\Phi \left(\frac{{\gamma }_{\lambda }}{{\delta }_{\lambda }}\right)}\varphi \left(\frac{{\gamma }_{\lambda }}{{\delta }_{\lambda }}\right)\frac{{\gamma }_{\lambda }}{\sqrt{2\pi }}\exp \left(-\frac{{\gamma }_{\lambda }^2}{2{\delta }_{\lambda }^2}\right)+{\gamma }_{\lambda }=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\lambda }}_i\hfill \\ 1-\frac{1}{\Phi \left(\frac{{\gamma }_{\lambda }}{{\delta }_{\lambda }}\right)}\varphi \left(\frac{{\gamma }_{\lambda }}{{\delta }_{\lambda }}\right)\frac{{\gamma }_{\lambda }^2}{\sqrt{2\pi }{\delta }_{\lambda }^2}\exp \left(-\frac{{\gamma }_{\lambda }^2}{2{\delta }_{\lambda }^2}\right)-{\gamma }_{\lambda }^2=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\lambda }}_i^2-\frac{2{\gamma }_{\lambda }}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\lambda }}_i\hfill \end{array}$$

$$\{\begin{array}{c}\frac{1}{\Phi \left(\frac{{\gamma }_{\sigma }}{{\delta }_{\sigma }}\right)}\varphi \left(\frac{{\gamma }_{\sigma }}{{\delta }_{\sigma }}\right)\frac{{\gamma }_{\sigma }}{\sqrt{2\pi }}\exp \left(-\frac{{\gamma }_{\sigma }^2}{2{\delta }_{\sigma }^2}\right)+{\gamma }_{\sigma }=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\sigma }}_i\hfill \\ 1-\frac{1}{\Phi \left(\frac{{\gamma }_{\sigma }}{{\delta }_{\sigma }}\right)}\varphi \left(\frac{{\gamma }_{\sigma }}{{\delta }_{\sigma }}\right)\frac{{\gamma }_{\sigma }^2}{\sqrt{2\pi }{\delta }_{\sigma }^2}\exp \left(-\frac{{\gamma }_{\sigma }^2}{2{\delta }_{\sigma }^2}\right)-{\gamma }_{\sigma }^2=\frac{1}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\sigma }}_i^2-\frac{2{\gamma }_{\sigma }}{n_{\mathrm{d}}}\sum _{i=1}^{n_{\mathrm{d}}}{\widehat{\sigma }}_i\hfill \end{array}$$

where $\widehat{\alpha }=\left({\widehat{\alpha }}_1,{\widehat{\alpha }}_2,...,{\widehat{\alpha }}_{n_{\mathrm{d}}}\right)$, $\widehat{\mu }=\left({\widehat{\mu }}_1,{\widehat{\mu }}_2,...,{\widehat{\mu }}_{n_{\mathrm{d}}}\right)$, $\widehat{\lambda }=\left({\widehat{\lambda }}_1,{\widehat{\lambda }}_2,...,{\widehat{\lambda }}_{n_{\mathrm{d}}}\right)$ and $\widehat{\sigma }=\left({\widehat{\sigma }}_1,{\widehat{\sigma }}_2,...,{\widehat{\sigma }}_{n_{\mathrm{d}}}\right)$ mean the training results of the corresponding model parameters.

Bayesian inference method is applied to estimate the model parameters. Generally, the parameters’ priors can be set based on the engineering experience and related information when applying Bayesian inference method. Furthermore, the non-informative prior also can be adopted for the conditions lacking related information. For the simulation and case studies presented in this paper, the priors of the hyper parameters are set based on TN distribution and the above MLEs. Furthermore, similarly as Ref. [18], the variances are set to be one-tenth of the MLEs of the hyper parameters for the simulation and case studies, whereas in other cases, it can be set based on engineering experience. The prior distributions of the hyper parameters for the simulation study in Section 4 and case study in Section 5, π(θ^H), are set as Eq. (24).

$$\{\begin{array}{c}{\gamma }_{\alpha }\sim TN\left({\widehat{\gamma }}_{\alpha },0.01{\widehat{\gamma }}_{\alpha }^2,0,+\infty \right),{\delta }_{\alpha }\sim TN\left({\widehat{\delta }}_{\alpha },0.01{\widehat{\delta }}_{\alpha }^2,0,+\infty \right)\hfill \\ {\gamma }_{\mu }\sim TN\left({\widehat{\gamma }}_{\mu },0.01{\widehat{\gamma }}_{\mu }^2,0,+\infty \right),{\delta }_{\mu }\sim TN\left({\widehat{\delta }}_{\mu },0.01{\widehat{\delta }}_{\mu }^2,0,+\infty \right)\hfill \\ {\gamma }_{\lambda }\sim TN\left({\widehat{\gamma }}_{\lambda },0.01{\widehat{\gamma }}_{\lambda }^2,0,+\infty \right),{\delta }_{\lambda }\sim TN\left({\widehat{\delta }}_{\lambda },0.01{\widehat{\delta }}_{\lambda }^2,0,+\infty \right)\hfill \\ {\gamma }_{\sigma }\sim TN\left({\widehat{\gamma }}_{\sigma },0.01{\widehat{\gamma }}_{\sigma }^2,0,+\infty \right),{\delta }_{\sigma }\sim TN\left({\widehat{\delta }}_{\sigma },0.01{\widehat{\delta }}_{\sigma }^2,0,+\infty \right)\hfill \end{array}$$

where the MLEs of the hyper parameters ${\widehat{\theta }}^{\mathrm{H}}=\left({\widehat{\gamma }}_{\alpha },{\widehat{\delta }}_{\alpha },{\widehat{\gamma }}_{\mu },{\widehat{\delta }}_{\mu },{\widehat{\gamma }}_{\lambda },{\widehat{\delta }}_{\lambda },{\widehat{\gamma }}_{\sigma },{\widehat{\delta }}_{\sigma }\right)$.

Based on the above hyper parameters’ priors, the Bayesian inference for the hyper parameters of the testing sample population, Eq. (25), is performed to estimate the hyper parameters for the testing sample population by solving Eq. (26) with Markov Chain Monte Carlo (MCMC) method.

$$p\left({\theta }^{\mathrm{H}}|\widetilde{ANN},Z_{1:n_{\mathrm{d}}},T_{1:n_{\mathrm{d}}}\right)=\frac{L\left(Z_{1:n_{\mathrm{d}}}|\widetilde{ANN},T_{1:n_{\mathrm{d}}},{\theta }^{\mathrm{H}}\right)\pi \left({\theta }^{\mathrm{H}}\right)}{\underset{{\theta }^{\mathrm{H}}}{\int }L\left(Z_{1:n_{\mathrm{d}}}|\widetilde{ANN},T_{1:n_{\mathrm{d}}},{\theta }^{\mathrm{H}}\right)\pi \left({\theta }^{\mathrm{H}}\right)d{\theta }^{\mathrm{H}}}$$

$$p\left({\theta }^{\mathrm{H}}|\widetilde{ANN},Z_{1:n_{\mathrm{d}}},T_{1:n_{\mathrm{d}}}\right)\propto L_{\mathrm{p}}\left(Z_{1:n_{\mathrm{d}}}|\widetilde{ANN},T_{1:n_{\mathrm{d}}},{\theta }^{\mathrm{H}}\right)\pi \left({\theta }^{\mathrm{H}}\right)$$

where $p\left({\theta }^{\mathrm{H}}|\widetilde{ANN},Z_{1:n_{\mathrm{d}}},T_{1:n_{\mathrm{d}}}\right)$ is the hyper parameters’ posterior distributions for the testing sample population. $L\left(Z_{1:n_{\mathrm{d}}}|\widetilde{ANN},T_{1:n_{\mathrm{d}}},{\theta }^{\mathrm{H}}\right)$ is the likelihood function for the testing sample population, $Z_{1:n_{\mathrm{d}}}$ and $T_{1:n_{\mathrm{d}}}$, and given by

$$L\left(Z_{1:n_{\mathrm{d}}}|\widetilde{ANN},T_{1:n_{\mathrm{d}}},{\theta }^{\mathrm{H}}\right)=\prod _{i=1}^{n_{\mathrm{d}}}\underset{{\theta }_i}{\int }f\left(Z_i|{\theta }_i,\widetilde{ANN},T_i\right)f_{\text{TN}}\left({\theta }_i|{\theta }^{\mathrm{H}}\right)d{\theta }_i$$

where f_TN(θ_i|θ^H) means the PDF of the parameter vector and given by

$$f_{\text{TN}}\left({\theta }_i|{\theta }^{\mathrm{H}}\right)=f_{\text{TN}}\left({\alpha }_i|{\gamma }_{\alpha },{\delta }_{\alpha }^2,0,+\infty \right)f_{\text{TN}}\left({\mu }_i|{\gamma }_{\mu },{\delta }_{\mu }^2,0,+\infty \right)f_{\text{TN}}\left({\lambda }_i|{\gamma }_{\lambda },{\delta }_{\lambda }^2,0,+\infty \right)f_{\text{TN}}\left({\sigma }_i|{\gamma }_{\sigma },{\delta }_{\sigma }^2,0,+\infty \right)$$

Based on Bayesian inference and MCMC method, the hyper parameters for the testing sample population are estimated and indicated as ${\overline{\theta }}^{\mathrm{H}}$. Given the failure threshold D, the reliability function and the PDF of the lifetime distribution for the testing sample population are given by Eqs. (29) and (30) [46].

$$R_{\mathrm{P}}\left(t|{\overline{\theta }}^{\mathrm{H}},\widetilde{ANN},D\right)=\underset{\theta }{\int }\Phi \left(\frac{D-\mu \widetilde{ANN}\left(\alpha t\right)}{\lambda \sqrt{\mu \widetilde{ANN}\left(\alpha t\right)}}\right)-\exp \left(\frac{2D}{{\lambda }^2}\right)\Phi \left(-\frac{D+\mu \widetilde{ANN}\left(\alpha t\right)}{\lambda \sqrt{\mu \widetilde{ANN}\left(\alpha t\right)}}\right)f_{\text{TN}}\left(\theta |{\overline{\theta }}^{\mathrm{H}}\right)d\theta$$

$$f_{\text{PT}}\left(t|{\overline{\theta }}^{\mathrm{H}},\widetilde{ANN},D\right)=\underset{\theta }{\int }\sqrt{\frac{D^2}{2\pi {\lambda }^2{\left(\mu \widetilde{ANN}\left(\alpha t\right)\right)}^3}}\exp \left[-\frac{{\left(\mu \widetilde{ANN}\left(\alpha t\right)-D\right)}^2}{2{\lambda }^2\mu \widetilde{ANN}\left(\alpha t\right)}\right]\frac{d\mu \widetilde{ANN}\left(\alpha t\right)}{dt}{f}_{\text{TN}}\left(\theta |{\overline{\theta }}^{\mathrm{H}}\right)d\theta$$

Besides the reliability and lifetime of the testing sample population, the individual reliability and lifetime are also noteworthy in engineering practice. When evaluating the model parameters for the monitoring individual, the model parameters’ priors, π(θ), can be set based on the above estimated hyper parameters, as

$$\{\begin{array}{c}\alpha \sim TN\left({\overline{\gamma }}_{\alpha },{\overline{\delta }}_{\alpha }^2,0,+\infty \right),\mu \sim TN\left({\overline{\gamma }}_{\mu },{\overline{\delta }}_{\mu }^2,0,+\infty \right)\hfill \\ \lambda \sim TN\left({\overline{\gamma }}_{\lambda },{\overline{\delta }}_{\lambda }^2,0,+\infty \right),\sigma \sim TN\left({\overline{\gamma }}_{\sigma },{\overline{\delta }}_{\sigma }^2,0,+\infty \right)\hfill \end{array}$$

Based on the above model parameters’ priors, the posterior distributions of the model parameters for the monitoring individual, Z and T, are estimated by Bayesian inference method, Eq. (32). The posterior distributions can be obtained based on Eq. (33) by MCMC method.

$$p\left(\theta |\widetilde{ANN},Z,T\right)=\frac{L_{\mathrm{I}}\left(Z|\theta ,\widetilde{ANN},T\right)\pi \left(\theta \right)}{\underset{\theta }{\int }{L}_{\mathrm{I}}\left(Z|\theta ,\widetilde{ANN},T\right)\pi \left(\theta \right)d\theta }$$

$$p\left(\theta |\widetilde{ANN},Z,T\right)\propto L_{\mathrm{I}}\left(Z|\theta ,\widetilde{ANN},T\right)\pi \left(\theta \right)$$

where $p\left(\theta |\widetilde{ANN},Z,T\right)$ means the posterior distributions of the model parameters for the monitoring individual. $L_{\mathrm{I}}\left(Z|\theta ,\widetilde{ANN},T\right)$ is the likelihood function for the monitoring individual and given by

$$L_{\mathrm{I}}\left(Z|\theta ,\widetilde{ANN},T\right)=f\left(Z|\theta ,\widetilde{ANN},T\right)$$

Based on the estimation results of the model parameters, $\overline{\theta }=\left(\overline{\alpha },\overline{\mu },\overline{\lambda },\overline{\sigma }\right)$, the reliability function and the PDF of the lifetime distribution for the monitoring individual are given by Eq. (35) and Eq. (36), respectively.

$$R_{\mathrm{I}}\left(t|\overline{\lambda },\overline{\Lambda }\right)=\Phi \left(\frac{D-\overline{\Lambda }\left(t\right)}{\overline{\lambda }\sqrt{\overline{\Lambda }\left(t\right)}}\right)-\exp \left(\frac{2D}{{\overline{\lambda }}^2}\right)\Phi \left(-\frac{D+\overline{\Lambda }\left(t\right)}{\overline{\lambda }\sqrt{\overline{\Lambda }\left(t\right)}}\right)$$

$$f_{\text{IT}}\left(t|\overline{\lambda },\overline{\Lambda },D\right)=\sqrt{\frac{D^2}{2\pi {\overline{\lambda }}^2{\left(\overline{\Lambda }\left(t\right)\right)}^3}}\exp \left[-\frac{{\left(\overline{\Lambda }\left(t\right)-D\right)}^2}{2{\overline{\lambda }}^2\overline{\Lambda }\left(t\right)}\right]\frac{d\overline{\Lambda }\left(t\right)}{dt}$$

where $\overline{\Lambda }\left(t\right)$ means the estimated mean function. Namely, $\overline{\Lambda }\left(t\right)=\overline{\mu }\widetilde{ANN}\left(\overline{\alpha }t\right)$.

The presented reliability estimation method can be applied to estimate the performances of the testing sample population and individual, including lifetime and reliability, considering the measurement error and individual difference among the units. Fig. 2 shows the procedure of the proposed method. Three steps need to be performed. First, the ANN supported mean function and the model parameters are obtained based on the training method shown in Section 3.1. Second, the training results of the model parameters are used to set the hyper parameters’ priors by MLE based method, Section 3.2. Finally, based on the above priors, the hyper parameters for the testing sample population are evaluated by Bayesian inference and used to evaluate the reliability and lifetime of the testing sample population, Section 3.3. The model parameters for the monitoring individual are estimated by Bayesian inference and used to infer the individual reliability and lifetime, in which the model parameters’ priors are set based on the above estimated hyper parameters. The corresponding algorithm is shown in Algorithm 1.

Algorithm 1. Algorithm of the proposed method.