Reliability Analysis of Small-Sample Failure Data for Random Truncation High-Voltage Relay

Abstract

In order to model and evaluate the reliability of long-life high-voltage relays with small-sample fault data characteristics, a reliability analysis method integrating average rank, the minimum mean square distance empirical distribution function, and total least squares estimation is proposed. In the random truncation experiment, considering the influence of random truncation data, the average rank method is used to correct the rank of small-sample fault data; then, the optimal empirical distribution function for small-sample fault data is obtained through the minimum average square distance, which can overcome the impact of small-sample fault data randomness. Under the assumption of the Weibull distribution model, the total least squares estimation method is used for reliability model parameter estimations, and the linear correlation coefficient and d-test method are used for model hypothesis testing. If two or more distribution models pass the linear correlation coefficient test and the d-test simultaneously, the root mean square error and relative root mean square error are applied to determine the optimal reliability model. The effectiveness of this method is verified by comparing it with the maximum likelihood estimation method.

1. Introduction

High-voltage relays are key electrical components of new energy vehicle power batteries, and their reliability directly affects the maintenance of power batteries and vehicle functions [1]. Reliability modeling and evaluations are the basis for evaluating the level of product reliability, facilitating reliability screening and reliability design. Therefore, it is imperative to conduct research on the reliability evaluation of electrical components in high-voltage systems. In recent years, with the improvement of product research and development level and the growth of the reliability level of high-voltage relays, the fault occurrence of high-voltage relays has been decreasing. How to accurately model and evaluate the reliability under small-sample fault data is one of the hotspots in the reliability research of high-voltage relays currently.

According to a literature search, there are various methods to solve the problem of insufficient fault sample data, such as the sample expansion method [2], the Bayesian method [3], the machine learning method [4], etc. The sample expansion method mainly includes the bootstrap method and the virtual augmented sample method [5]. The bootstrap method is simple and effective. It is applicable to small-sample data of any distribution type and has a good expansion effect on small-sample data with a sample size greater than 10. The virtual augmented sample method is suitable for extremely small-sample data with a sample size of less than 10. Before sample expansion takes place, this method usually uses prior information of similar products and combines prior information with small-sample data for sample expansion. However, for new products with fast updates, it is difficult to have sufficient and accurate empirical data to utilize. The Bayesian estimation method uses historical information or prior information of similar products to estimate the parameters of the reliability model and corrects the parameters to obtain the posterior distribution of the reliability model parameters through the obtained small-sample data, thereby conducting a reliability evaluation [6,7]. Mc Neish [8] pointed out that the selection of Bayesian prior distribution in small-sample situations has a very sensitive impact on the evaluation results. Machine learning methods such as neural networks [9] and support vector machines (SVMs) [10] are also used to analyze reliability under a small-sample-data situation. However, neural networks have problems such as a slow convergence speed, inconsistent network structure selection, and local optimal solutions. SVMs are extremely sensitive to sample errors and lack adaptive ability, which can greatly reduce the robustness of reliability analyses.

The parameter estimation methods for reliability distribution models can be divided into the moment method, the maximum likelihood method, the differential evolution method, the Bayesian method, the grey model method, and the least squares estimation method [11,12]. Most of these methods do not have analytical solutions and require solving nonlinear equations or the maximum points of multivariate functions. The grey model method can estimate the three parameters of three-parameter Weibull distribution at once without iteration, and it avoids the problem of correlations existing between parameters and has a fast calculation speed, making it suitable for small-sample Weibull distribution data analysis [13]. However, the principle of this method limits the order of its fault data to integers, so it is not suitable for small-sample situations with missing data. The least squares estimation method is widely used due to its simple principle and analytical solution. It also has the advantage of being independent of the initial values. The assumption distribution model in engineering practice is usually a two-parameter Weibull distribution [14]. The principle of the least squares method for two-parameter Weibull distribution parameter estimation is to convert the Weibull distribution function into a linear regression model and obtain parameter estimates by linearly fitting the empirical distribution function values with fault data. It can be seen that different fault ranks and the selection of empirical distribution functions result in different parameter estimation results. Therefore, accurate fault ranks and finding the optimal empirical distribution function are crucial for the parameter estimation results of the least squares method [15,16].

This article focuses on the reliability modeling and evaluating of long-life high-voltage relays with small-sample fault data characteristics and proposes a reliability analysis method integrating average rank, the minimum mean square distance empirical distribution function, and total least squares estimation. In the random truncation experiment, considering the influence of random truncation data, the average rank method is used to correct the rank of small-sample fault data; then, the optimal empirical distribution function for small-sample fault data is obtained through the minimum average square distance. Under the assumption of the Weibull distribution model, the total least squares estimation method is used for reliability model parameter estimations, and the linear correlation coefficient and d-test method are used for model hypothesis testing. If two or more distribution models pass the linear correlation coefficient test and the d-test simultaneously, the root mean square error and relative root mean square error are applied to deter-mine the optimal reliability model. The effectiveness of this method is verified by comparing it with the maximum likelihood estimation method.

2. Principle of Random Truncation Small-Sample Reliability Analysis

In order to save time and cost, many products use a random truncation test method for reliability testing. The characteristics of fault data under a random truncation test can be represented in Figure 1.

Regarding highly reliable and long-life products, most products do not fail at the end of the experiment, resulting in a large amount of truncated data and a small amount of fault data. The existence of truncated data directly affects the fault rank. Therefore, fault rank corrections are required. Meanwhile, considering the impact of empirical distribution functions on parameter estimations, the minimum average square distance is used to select the optimal empirical distribution function. The average square distance does not depend on the value of the fault samples and the optimal empirical distribution function can be selected from the perspective of the overall average. Under the assumption of the Weibull distribution model, the total least squares estimation method is introduced for parameter estimation. Then, linear correlation coefficient and d-test method are used for model hypothesis testing. If two or more distribution models pass the linear correlation coefficient test and the d-test simultaneously, the root mean square error and relative root mean square error are applied to determine the optimal reliability model. A flow chart of the reliability analysis of random truncation small-sample fault data is shown in Figure 2.

2.1. Basic Principle of Average Rank Method

Assuming five samples of the same type are subjected to random truncation tests, the failure time and truncation time are arranged in ascending order, denoted as $t_1\le t_2+\le t_3\le t_4\le t_5$, where “$+$” represents right truncated data. If the test continues until the No.2 product fails, then the fault time $t_2$ may occur at any time within the interval $\left(t_2+,t_3\right),\left(t_3,t_4\right)\left(t_4,t_5\right)\left(t_5,\infty \right)$. It can be considered that the failure probability of the No. 2 product is the same within these four time intervals. The fault rank of $t_1$ is 1 regardless of which time interval $t_2$ occurred in. The fault rank of $t_3$ is between 2 and 3, with a probability 3/4 of 2 and a probability 1/4 of 3. From the perspective of equal probability, the fault rank of $t_3$ is $2\times 3/4+3\times 1/4=9/4$. Similarly, the fault rank of $t_4$ and $t_5$ are $3\times 2/4+4\times 2/4=14/4$ and $4\times 1/4+5\times 3/4=19/4$, respectively.

The above fault rank calculation process can be generalized as follows. Assuming $n$ samples are subjected to random truncation experiments, of which $p$ samples are faulted and $n-p$ samples are right truncated, we arrange all the $n$ data in ascending order, and their ranks are numbered from 1 to $n$, denoted as $j$. Then, arrange $p$ fault data in ascending order, and the fault rank of the $i$th fault data is calculated using Equation (1):

$$r_i=r_{i-1}+\left(n+1-r_{i-1}\right)/\left(n+2-j\right)$$

(1)

where $r_i$ is the fault rank of the $i$th fault data, $1\le i\le p$, $r_0=0$.

2.2. Optimal Selection of Empirical Distribution Functions

Assuming $X_1,X_2,\cdots ,X_n$ is an independent and identically distributed sample from population $X$, the $i$th order statistic of this sample is denoted as $X_{\left(i\right)}$. $x_1,x_2,\cdots ,x_n$ are observations of $X_1,X_2,\cdots ,X_n$. We arrange them in ascending order, $x_{\left(1\right)}\le x_{\left(2\right)}\le \cdots \le x_{\left(n\right)}$, where $x_{\left(i\right)}$ is $i$th order statistic. The commonly defined empirical distribution function $Q_n\left(x\right)$ is as follows:

$$Q_n\left(x\right)=\left\{\begin{array}{c}0,x\le x_{\left(1\right)};\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \frac{k}{n},x_{\left(k\right)}\le x\le x_{\left(k+1\right)},k=1,2,\cdots ,n-1;\\ 1,x_{\left(n\right)}\le x.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(2)

where $n$ is the number of the samples; $k$ is the ordinal number of order statistics.

From Equation (2), it can be seen that $Q_n\left(x\right)$ is a non-decreasing left continuous function with $Q_n\left(-\infty \right)=0,Q_n\left(+\infty \right)=1$. Another empirical distribution function, $G_n\left(x\right)$, recommended by reference [17] is as follows:

$$G_n\left(x\right)=\left\{\begin{array}{l}0,\hspace{1em}x\le x_{\left(1\right)};\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ \frac{k-0.375}{n+0.25},\hspace{1em}{x}_{\left(k\right)}\le x\le x_{\left(k+1\right)},k=1,2,\cdots ,n-1;\\ \frac{n-0.375}{n+0.25},\hspace{1em}{x}_{\left(n\right)}\le x.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(3)

According to (2) and (3), the general expression for the empirical distribution function is as follows:

$$M_n\left(x\right)=\left\{\begin{array}{c}{m}_0,\hspace{1em}x\le x_{\left(1\right)};\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ m_k,\hspace{1em}{x}_{\left(k\right)}\le x\le x_{\left(k+1\right)},k=1,2,\cdots ,n-1;\\ m_n,\hspace{1em}{x}_{\left(n\right)}\le x.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(4)

In the study of probability graphs, there are many expressions for $m_k$, and the commonly used ones of $m_k$ include expectation estimation $m_k=k/\left(n+1\right)$, midpoint estimation $m_k=\left(k-0.5\right)/n$, median rank estimation $m_k=\left(k-0.3\right)/\left(n+0.4\right)$, and approximate median rank estimation $m_k=\left(i+1\right)/\left(n+2\right)$, where $m_0=0$ [18].

The selection of $m_k$ depends on the approximation degree between the empirical distribution $M_n\left(x\right)$ and the overall distribution $F\left(x\right)$. The maximum linear distance $d_n\left(M\right)=\sup \left|M_n\left(x\right)-F\left(x\right)\right|$ is intuitive and has practical significance. However, it only considers the points of extreme values and does not examine the approximation degree between $M_n\left(x\right)$ and $F\left(x\right)$ from an overall perspective. The square distance $L\left(M_n,F\right)={{\displaystyle \int \left(M_n-F\right)}}^{2}dF$ examines the approximation degree from an overall perspective and reflects the difference between $M_n\left(x\right)$ and $F\left(x\right)$ more comprehensively. However, square distance depends on the value of the sample, resulting in randomness in the outcome. Therefore, in the average sense, this article introduces the average square distance to measure the overall approximation degree of the sample, which can overcome the impact of small-sample fault data randomness [18]. The average square distance is calculated as follows:

$$R\left(M_n\right)=E\left(L\left(M_n,F\right)\right)={\displaystyle \int \dots {\displaystyle \int L\left(M_n,F\right)d{P}_n\left(x_1,x_2,\cdots ,x_n\right)}}$$

(5)

where $P_n\left(x_1,x_2,\cdots ,x_n\right)$ is the joint distribution function of $\left(x_{\left(1\right)},x_{\left(2\right)},\cdots ,x_{\left(n\right)}\right)$.

Therefore, the $m_k$ with the minimum average square distance is the optimal empirical distribution function.

2.3. Total Least Squares Parameter Estimation

The principle of least squares parameter estimation is to obtain the parameter estimation value of the linear model $y=a_{1}x+b_1$ by minimizing the sum of squared residuals, i.e., $min{\displaystyle \sum _{i=1}^p{\left[y_i-\left(a_1{x}_i+b_1\right)\right]}^2}$. At this point, the estimated values of the regression coefficients ${\widehat{a}}_1$ and ${\widehat{b}}_1$ are calculated according to Equation (6) [19].

$$\left\{\begin{array}{l}{\widehat{a}}_1=\frac{{\displaystyle \sum _{i=1}^p\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{{\displaystyle \sum _{i=1}^p{\left(x_i-\overline{x}\right)}^2}}\\ {\widehat{b}}_1=\overline{y}-{\widehat{a}}_1\overline{x}\end{array}\right.$$

(6)

where $\overline{x}={\displaystyle \sum _{i=1}^p{x}_i}/p,\overline{y}={\displaystyle \sum _{i=1}^p{y}_i}/p$, and $p$ is the number of sample points used for fitting.

The total least squares method [20] is used to obtain the parameter estimation value of the linear model $y=a_{2}x+b_2$ by minimizing the sum of squared vertical residuals, i.e., $min\frac{1}{1+a^2}{\displaystyle \sum _{i=1}^p{\left[y_i-\left(a_2{x}_i+b_2\right)\right]}^2}$. The estimated values of the regression coefficients ${\widehat{a}}_2$ and ${\widehat{b}}_2$ are calculated according to Equation (7).

$$\left\{\begin{array}{l}{\widehat{a}}_2=-c\pm \sqrt{1+c^2}\\ {\widehat{b}}_2=\overline{y}-{\widehat{a}}_2\overline{x}\end{array}\right.$$

(7)

where $c=\frac{{\displaystyle \sum _{i=1}^p{\left(x_i-\overline{x}\right)}^2-}{\displaystyle \sum _{i=1}^p{\left(y_i-\overline{y}\right)}^2}}{2{\displaystyle \sum _{i=1}^p\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}$; “$\pm$” in ${\widehat{a}}_2$ depends on the distribution of $\left(x_i,y_i\right)\left(i=1,2,\cdots ,p\right)$. If the sample covariance is greater than 0, we take “$+$”; if it is less than 0, we take “$+$“.

The common reliability distribution function in engineering is the Weibull distribution, $F\left(t\right)=1-\exp \left[-{\left(t/\eta \right)}^{\beta }\right]$, where $x=\ln t$, $y=\ln \left[\ln \left(1/\left(1-F\left(t\right)\right)\right)\right]$ can be transformed into a linear model, $y=ax+b$ [21]. In terms of practical significance, it can be seen that $x_i$ and $y_i$ are positively correlated, and the covariance is greater than 0. Therefore, in Formula (7), ${\widehat{a}}_2=-c+\sqrt{1+c^2}$. At this point, the Weibull distribution parameter estimations are as follows.

$$\left\{\begin{array}{l}\widehat{\eta }=\exp \left(-{\widehat{b}}_2/{\widehat{a}}_2\right)\\ \stackrel{⌢}{\beta }={\widehat{a}}_2=-c+\sqrt{1+c^2}\end{array}\right.$$

(8)

where $\stackrel{⌢}{\beta }$ is the Weibull distribution shape parameter, and $\widehat{\eta }$ is the Weibull distribution scale parameter.

2.4. Model Hypothesis Testing

2.4.1. Linear Correlation Test

Using linear correlation coefficient $\widehat{\rho }$ to test the linear correlation between linear models, the formula for calculating the correlation coefficient is as follows [22]:

$$\widehat{\rho }=\frac{{\displaystyle \sum _{i=1}^p{x}_i{y}_i-n\overline{x}\hspace{0.17em}\overline{y}}}{\sqrt{\left({\displaystyle \sum _{i=1}^p{x}_i^2-p{\overline{x}}^2}\right)\left({\displaystyle \sum _{i=1}^p{y}_i^2-p{\overline{y}}^2}\right)}}$$

(9)

where $x_i=\ln t_i$, $y_i=\ln \left[\ln \left(1/\left(1-F\left(t_i\right)\right)\right)\right]$, $t_i$ is time of failure, $\overline{x}$ and $\overline{y}$ are average values of $x_i,y_i$ respectively, and $p$ is the number of sample points.

When $\widehat{\rho }>{\rho }_a$, it was believed that there was a linear correlation between $x$ and $y$, indicating the validity of the parameter estimation. ${\rho }_a$ is the minimum value of the correlation coefficient ($\alpha$ is the significance level, and the value of the significance level in engineering is usually 0.1, that is, $\alpha =0.1$) and ${\rho }_a$ is obtained by the critical values table of linear correlation coefficients.

2.4.2. D-Test

To test the Weibull distribution hypothesis and the fitting effect of the distribution, the d-test is required. The d-test, also known as the Kolmogorov test, arranges $p$ fault data in ascending order, calculates the distribution function $\widehat{F}(t_i)$ corresponding to each data $t_i$ based on parameter estimation, and compares $\widehat{F}(t_i)$ with the empirical distribution function $M\left(t_i\right)$. The d-test statistic is as follows:

$$D_p=\underset{-\infty <x<+\infty }{\sup }\left|\stackrel{\wedge }{F}\left(t_i\right)-M\left(t_i\right)\right|=\max \left\{d_i\right\}$$

(10)

We compare $D_p$ with the critical value $D_{p,\alpha }$ and accept the original hypothesis if $D_p<D_{p,\alpha }$. Otherwise, we reject the original hypothesis [23]. $D_{p,\alpha }$ can be obtained via the d-test statistical table.

2.5. Model Selection and Reliability Evaluation

2.5.1. Model Selection

If two or more distribution models pass the linear correlation coefficient test and the d-test simultaneously, the root mean square error (RMSE) and relative root mean square error (MRMSE) are applied to determine the optimal reliability model, as shown in Formulas (11) and (12) [24].

$$RMSE=\sqrt{\frac{1}{p}{\displaystyle \sum _{i=1}^p{\left[\stackrel{\wedge }{F}\left(t_i\right)-M\left(t_i\right)\right]}^2}}$$

(11)

$$MRMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^p{\left[\stackrel{\wedge }{F}\left(t_i\right)-M\left(t_i\right)\right]}^2}}{{\displaystyle \sum _{i=1}^p{\left[M\left(t_i\right)\right]}^2}}}$$

(12)

where $\stackrel{\wedge }{F\left(t_i\right)}$ is the distribution function value of the distribution model obtained from the parameter estimation, and $M\left(t_i\right)$ is the empirical distribution function.

2.5.2. Reliability Assessment

After the optimal distribution model is determined, MTBF (mean time between failure) is used to evaluate product’s reliability [25].

The point estimation of MTBF is calculated using Equation (13):

$$MTBF={\displaystyle {\int }_0^{\infty }t}f(t)dt={\displaystyle {\int }_0^{\infty }t}dF(t)=tF(t){|}_0^{\infty }-{\displaystyle {\int }_0^{\infty }F(t)dt}$$

(13)

where $f\left(t\right)$ is the probability density function, and $F\left(t\right)$ is the failure distribution function.

3. Application Cases

Moreover, 56 electric vehicles with Type I high-voltage relays were subjected to random truncation reliability tests simultaneously. If the high-voltage relay failed, a new high-voltage relay was used as a replacement. A total of 22 electric vehicle samples generated eight relays fault data as shown in

Table 1. Fault data of relay reliability test (unit: km).

Relay Number	1–4	5–8	9–14	15–18	19–22
Fault occurrence time	50	100	500	2000	500
Total test time T	7920	11,700	4000	26,000	16,000
Number of faults	1	2	1	3	1

, with $T$ being the test time for each relay. The testing time for another 34 samples is shown in

Table 2. Reliability test truncation time (unit: km), with “+” representing truncated data.

Relay Number	23–26	27–30	31–34	35–40	41–44	45–48	49~52	53~56
Total test time T	5000+	6000+	7000+	2000+	3000+	30,000+	1000+	500+

.

According to Equation (1), the fault rank i of eight fault data can be calculated, as shown in

Table 3. Fault rank of relay test fault data (unit: km), with “+” representing truncated data.

$\mathbf{Data} {\mathit{t}}_{\mathit{j}}$	Order Number j (Ascending Order)	Fault Rank i (Equation (1))
50	1	1
100	3	3
500+	7
500	9	5.1379
1000+	13
2000+	19
2000	22	9.0419
>2000	64

.

According to Equation (4) and the five commonly used expressions of $m_k$, five different empirical distribution functions can be obtained, as shown in

Table 4. Empirical distribution function, with “+” representing truncated data.

$\mathbf{Data} {\mathit{t}}_{\mathit{j}}$	Fault Rank i	${\mathit{F}}_1\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}+1}{\mathit{n}+2}$	${\mathit{F}}_2\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.3}{\mathit{n}+0.4}$	${\mathit{F}}_3\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}}{\mathit{n}+1}$	${\mathit{F}}_4\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.5}{\mathit{n}}$	${\mathit{F}}_5\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.375}{\mathit{n}+0.25}$
50	1	0.030303	0.010870	0.015384	0.007812	0.009712
100	3	0.0606061	0.041925	0.046153	0.039062	0.040792
500+
500	5.1379	0.902998	0.075123	0.079044	0.072467	0.074015
1000+
2000+
2000	9.0419	0.15215	0.135744	0.139106	0.133467	0.134683
>2000

. $n=64$ with 8 fault data and 56 truncation data.

According to Equation (5), the average square distance of each empirical distribution function can be obtained, as shown in

Table 5. Average square distance.

Empirical Distribution Function	${\mathit{F}}_1\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}+1}{\mathit{n}+2}$	${\mathit{F}}_2\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.3}{\mathit{n}+0.4}$	${\mathit{F}}_3\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}}{\mathit{n}+1}$	${\mathit{F}}_4\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.5}{\mathit{n}}$	${\mathit{F}}_5\left({\mathit{t}}_{\mathit{i}}\right)=\frac{\mathit{i}-0.375}{\mathit{n}+0.25}$
$R\left(M_n\right)$	0.00260	0.00263	0.00252	0.00266	0.00252

.

From Table 5, it can be seen that the average square distance of $F_3\left(t_i\right)=\frac{i}{n+1}$ and $F_5\left(t_i\right)=\frac{i-0.375}{n+0.25}$ are almost the same and minimum.

Therefore, taking $F_3\left(t_i\right)=\frac{i}{n+1}$ and $F_5\left(t_i\right)=\frac{i-0.375}{n+0.25}$ as the optimal empirical distribution functions, assuming that the life distribution of high-voltage relays is a Weibull distribution, using the total least squares method, and according to Formulas (7) and (8), Weibull distribution parameters can be estimated. Then, according to Formulas (9) and (10), the linear correlation test and d-test are performed. The parameter estimation values and test values are shown in

Table 6. Weibull distribution parameter estimation of relay test fault data.

Empirical Distribution Function	β/η of Weibull Distribution	Linear Correlation Coefficient	D-Test Value	Root Mean Square Error/Relative Root Mean Square Error
$F_3\left(t_i\right)=\frac{i}{n+1}$	0.5347/70,623	0.952	0.0166	0.0116/0.1501
$F_5\left(t_i\right)=\frac{i-0.375}{n+0.25}$	0.6288/35,586	0.922	0.0293	0.0175/0.2361

.

When the significance level $\alpha =0.1$, the critical values of the linear correlation coefficient ${\widehat{\rho }}_{0.1}$ and the critical values of the d-test method $D_{p,\alpha }$ can be obtained through a table search, which are ${\widehat{\rho }}_{0.1}=0.9000$, $D_{p,\alpha }=0.51$, respectively. As shown in Table 6, both Weibull distribution functions obtained from the two empirical distribution functions pass the linear correlation coefficient test and the d-test. Therefore, the root mean square error and relative root mean square error methods are introduced to optimize the reliability model, the results of which are also shown in Table 6.

Table 6 shows that based on the root mean square error value and relative root mean square error, the priority order of the two empirical distribution functions is 3, 5. The conclusion drawn from the values of the root mean square error and the relative root mean square error is consistent; that is, the Weibull distribution model estimated based on $F_3\left(t_i\right)=\frac{i}{n+1}$ is the best. The optimal Weibull distribution parameters are as follows:

$$\left\{\begin{array}{l}\widehat{\eta }=\mathrm{70,623}\\ \stackrel{⌢}{\beta }=0.5347\end{array}\right.$$

The optimal Weibull distribution is

$$F\left(t\right)=1-\exp \left[-{\left(t/\mathrm{70,623}\right)}^{0.5347}\right]$$

According to Equation (13), the MTBF of the high-voltage relay is

$$\begin{array}{l}MTBF={\displaystyle {\int }_0^{\infty }t}f(t)dt={\displaystyle {\int }_0^{\infty }t}dF(t)=tF(t){|}_0^{\infty }-{\displaystyle {\int }_0^{\infty }F(t)dt}\\ =\eta \mathsf{\Gamma }\left(1+1/\beta \right)=\mathrm{125,741.04} \mathrm{km}\end{array}$$

Comparing the maximum likelihood estimation method with the method proposed in this paper, the maximum likelihood function of the Weibull distribution is:

$$\begin{array}{l}L\left(t_i;\eta ,\beta \right)={\displaystyle \prod _{i=1}^{r}f\left(t_i;\eta ,\beta \right)}{\displaystyle \prod _{i=1}^{N-r}R\left(t_i;\eta ,\beta \right)}\\ =\frac{{\beta }^r}{{\eta }^{r\beta }}{\left({\displaystyle \prod _{i=1}^r{t}^{\beta -1}}\right)}^{\beta -1}\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\displaystyle \sum _{i=1}^n{\left(t/\eta \right)}^{\beta }}\right]\end{array}$$

(14)

where $r$ is the number of failures and $N$ is the sum of the number of failures and truncations.

After logarithmic differentiation calculations, the estimated values of Weibull distribution parameters can be obtained as follows:

$$\left\{\begin{array}{l}\widehat{\eta }’=\mathrm{91,085}\\ \stackrel{⌢}{\beta }’=0.576\end{array}\right.$$

The $MTBF’=\eta ’\mathsf{\Gamma }\left(1+1/\beta ’\right)$ obtained by the maximum likelihood method is 144112.08 km.

The definition of average life in reference [26] is as follows:

$$\stackrel{\wedge }{\theta }=\frac{T*}{r}$$

(15)

where $\stackrel{\wedge }{\theta }$ is the point estimate of MTBF; $T*$ is the total working time of the products; and $r$ is the number of failures.

After calculation, $\stackrel{\wedge }{\theta }=\mathrm{61,560} \mathrm{km}$. We compared $MTBF’=\mathrm{144,112.08} \mathrm{km}$ obtained by the maximum likelihood estimation method and $MTBF=\mathrm{125,741.04} \mathrm{km}$ obtained by the method proposed in this paper with $\stackrel{\wedge }{\theta }$. It could be seen intuitively that the MTBF obtained by the method proposed in this paper was closer to $\stackrel{\wedge }{\theta }$, which demonstrated the effectiveness and superiority of the proposed method in small-sample reliability assessments.

4. Conclusions

(1) In order to solve the problem of inaccurate reliability modeling and evaluations caused by insufficient fault data of long-life high-voltage relays, this paper applies the average rank method to correct the influence of truncation data on the fault rank. Based on this, the minimum average square distance is used to select the optimal empirical distribution function, which avoids the randomness of small-sample fault data. After calculations, it is seen that average square distances of the empirical distribution function $F_3\left(t_i\right)=\frac{i}{n+1}$ and $F_5\left(t_i\right)=\frac{i-0.375}{n+0.25}$ are equal and minimum. Under the assumption of the Weibull distribution, the total least squares method is used to estimate Weibull distribution model parameters. The linear correlation coefficient test and d-test are used for model hypothesis testing and verifying the validity of the parameter estimations.

(2) In response to the optimal model selection problem for multiple reliability models that pass both the linear correlation test and d-test simultaneously, the root mean square error method and relative root mean square error are introduced, and the optimal Weibull distribution model are determined based on the principle of minimizing both.

(3) Based on the Weibull model, if the shape parameter is less than 1, it can be inferred that the product is in the early failure period. The reliability level of high-voltage relays is evaluated by estimating MTBF (mean time between failures). After carrying out calculations, it is found that compared with traditional methods, the MTBF obtained by the method proposed in this article is closer to the actual mean life of the product, indicating the effectiveness and superiority of the method proposed in this article. The method proposed in this article can lay the foundation for reliability research such as reliability design and screening, and can provide reference and guidance for the reliability modeling and life assessment of similar products.