Reliability Evaluation of Smart Meters Under Degradation-Shock Loads Based on Phase-Type Distributions

To address the problem that high-reliability and long-life complex products are susceptible to coupled fault mechanisms, and to address the challenge of constructing a macro fault behavior model for reliability evaluation, the reliability of a smart meter based on a Phase-Type (PH) distribution under a degradation-shock load is examined. This reliability assessment method considers the impacts of environmental stress and lightning shock on a smart meter. The performance degradation process caused by environmental stress is described by a Wiener model, which is converted into a representation of a PH distribution in a subsequent step. A lightning impulse can cause performance degradation or direct failure of a smart meter. The arrival time of the lightning, the magnitude of the lightning, and the damage caused by the lightning are analyzed separately, and then, the obtained model is expressed in the form of a PH distribution. Finally, a reliability evaluation model of a smart meter based on the PH model under a degradation-shock load is established. The accuracy of the model is verified by comparing the model results with Monte Carlo simulation results.


PH
Phase-Type HErD hyper-Erlang distribution Q Transition rate matrix of PH distribution a Initial probability of each state in a PH distribution S Transition rate between transient states S 0 Transition rate from the transient states to absorbing state f (x) The density function of X F(x) The distribution function of X α Initial probability of each Erlang branch The associate editor coordinating the review of this manuscript and approving it for publication was Cheng Qian .
λ Scale parameter of an Erlang distribution r Number of phases (states) of each Erlang branch P 1 , P 2 , P 3 the probability that the shock belongs to the damage, fatal and safety zone Threshold of safety zone D 1 Threshold of damage zone Y j The damage caused by the j-th shock W i The magnitude of the i-th shock Z k The time between two shocks X i the damage caused by environmental stress for one year. H The failure threshold of performance degradation N 1 (t) , N 2 (t) , N 3 (t) the number of shocks in the damage, fatal and safety zones prior to time t F D tm (H ) The probability of the damage caused by environmental stress and shock not exceeding the failure threshold at time t

I. INTRODUCTION
The smart meter is one of the most critical elements of smart grids. Its reliability directly affects the reliability and safety of the power supply for thousands of households. However, under the combined effect of temperature, humidity, electricity and other stressors, the internal components of the smart meter gradually deteriorate, leading to reduced reliability of the entire meter [1]. Therefore, reliability estimation is a significant technique to examine the quality of smart meters manufactured for power suppliers. The smart meter consists of different functional modules, including a power module, metering module, control module, display module, storage module, and communication module [2]. He et al. [3] used the failure mode mechanism and effect analysis (FMMEA) method to analyze the failure mechanism of a smart meter metering module and obtained the main sensitive stresses of the metering module: temperature, humidity, electrical stress and mechanical stress. However, most articles on the reliability of smart meters are based on experimental data. First, research on degradation (the performance of smart meters gradually degrades when they are affected by environmental stresses such as temperature and humidity) is as follows: Ye et al. [4] presented an experimental method to evaluate the smart meter acceleration life based on a degenerate trajectory model. Yang et al. [1] used the Weibull distribution to fit data obtained through an accelerated degradation test to obtain the reliability of smart meters. Xu et al. [5] developed a multivariate degradation modeling method via vine copulas to estimate the reliability of products with multiple PCs reflecting degradation states. Regarding research on shock (smart meters can also fail under external shocks such as lightning), according to the data records of fault smart meters in Chongqing in the second half of 2018, 334 of the 192,670 smart meters were found to fail due to lightning strikes. The impact of lightning on smart meters therefore cannot be ignored. Kazuyuki Ishimoto et al. [6] calculated the lightning damage rate of a smart meter due to a magnetic field by lightning surge analysis. In actual use, smart meters are subject to the combined effects of multiple stresses, so degradation-shock modeling needs to be studied in the reliability assessment of smart meters. Although the degradation-shock model has not been studied much in this context, it has long been studied in other fields. Lehman [7] considered a class of degradationthreshold-shock (DTS) models that provide a rich conceptual framework for studying degradation and shock. Based on this research, some researchers have noted a correlation between the degradation process and the shock process, and they have divided the impact into two types: The first type assumes that arriving shocks affect the degradation processes in some aspects, such as sudden increments [8], [9], shifting thresholds [10] and degradation rate variation [11], [12]. In addition, Jiang et al. [13] considered that the shocks can be divided into different zones, where each zone has a different effect on the degradation process. The second type assumes that the degradation process affects the shock processes in terms of the arrival intensity [14], [15] or destructive probability [16].
Phase-Type (PH) distributions are a class of mixture distributions that contain many well-known distributions as members, such as the exponential and Erlang distributions and their combinations. Using the maximum likelihood or expectation maximization (EM) algorithm, any distribution or data can be approximated by the PH distribution, which can reduce the difficulty of analysis and modeling by converting different types of distributions into PH representations [17]. Furthermore, matrix-geometric properties of these representations allow problems involving the calculation of convolutions to be handled [18]. PH distributions have been widely used in reliability analysis modeling. However, most of the literature focuses on modeling the shock process. Riascos-Ochoa et al. [19] proposed a reliability assessment model of a cumulative impact degradation system based on a PH distribution. Segovia and Labeau [20] concluded that when the shock amplitude does not reach the shock threshold, the values of the initial vector of the PH distribution change, increasing the probability that the system will enter a dangerous state at the next moment. Gong et al. [21] developed a generalized run shock model with two thresholds. Zhao et al. [22] proposed a multistate shock model in which three types of shocks were considered. In the shock model proposed by Ranjkesh et al. failure of the system means that the damage caused by the impact exceeds the threshold or that the interval between shock arrivals is less than the threshold [23]. However, in actual engineering, the effects of cumulative degradation and random shocks should be analyzed simultaneously. To date, few degradation-shock models based on PH distributions have been studied. Li et al. [24] discussed the modeling issue of a system subjected to a degradation-shock load based on the PH distribution, but they did not consider hard failure. There is a constraint relationship in which the shock magnitude is proportional to the shock damage.
In this paper, a reliability assessment method for smart meters based on a PH distribution under a degradation-shock load is proposed, and its framework is shown in Fig. 1. The degradation-shock model of this paper considers performance degradation caused by temperature and humidity as well as shock damage and direct failure caused by lightning. Using the PH fitting method, it is convenient to obtain various parameters of the model to evaluate the reliability of the smart meter. Finally, through comparison with the Monte Carlo method, the accuracy of the model constructed in this paper is verified.
The rest of the paper is organized as follows. Section II introduces the background of the PH distribution and its fitting method. Section III analyzes the arrival time, magnitude of lightning and damage caused by lightning. Section IV establishes a degradation-shock model based on the PH distribution under the combined action of environmental stresses and lightning shock. Section V discusses the conclusions.

II. PH DISTRIBUTION A. DEFINITION OF PH DISTRIBUTIONS
The PH distribution was introduced by Neuts [25], [26] and is created by a convolution or mixture of exponential distributions. Mathematically, it is defined as the continuous distribution of the time until absorption in a continuous time Markov chain with an absorbing state. To formally introduce this model, consider a continuous time Markov chain with m + 1 states, where states 1, . . . , m are transient states and state m + 1 is the absorbing state. Furthermore, let the process have an initial probability starting in any m + 1 phases given by the probability vector (a 0 , a), where a 0 is a scalar and a is a 1 × m vector.
This process can be written in the form of the following transition rate matrix: where 0 is a 1 ×m row vector consisting of zeros and 1 is an m × 1 column vector consisting of ones. Here, S is an m × m matrix, where m is commonly called the phase size. In this process, the distribution of the time X at which the absorption state is reached is referred to as a PH distribution and is expressed as PH (a, S). The density function and distribution function of X are given by where exp (·) is the matrix index. It is generally assumed that the probability of starting processing in the absorbing state is zero (i.e., a 0 = 0).
The following examples show how these distributions can be expressed as PH distributions.
Example 1 (Exponential Distribution): The simplest example of a PH distribution is an exponential distribution with a parameter of λ. The parameters of the PH distribution are S = −λ and a = 1.

Example 2 (Erlang Distribution):
The Erlang distribution has two parameters, namely, the shape parameter k > 0 and the rate parameter λ > 0, and is sometimes expressed as E (k, λ). The Erlang distribution can be written as a PH distribution by using a diagonal k-k matrix of diagonal elements -λ and superdiagonal elements λ, where the probability of starting in state 1 is equal to one.
PH distributions are dense in the set of continuous density functions with support on [0, ∞), which means that there is a PH distribution arbitrarily close to any continuous distribution. There are currently many effective algorithms [17], [27], [28] to fit the PH distribution to any (positive) data set.

2) CLOSURE UNDER CONVOLUTIONS
The convolution of PH distributions is also a PH distribution. In other words, the sum of two or more independent PH-distributed random variables is also PH distributed. If X 1 and X 2 are two independent random variables whose distribution is PH (τ 1 , T 1 ) of order n1 and PH (τ 2 , T 2 ) of order n 2 , then the sum is distributed as PH (a, S) of order n 1 + n 2 with where T 1 · 1 + t 1 = 0. Note that the sum of k independent PH-distributed random variables is also PH distributed, and its PH representation can be obtained by applying (4).

C. FITTING ALGORITHM OF PH DISTRIBUTIONS
The maximum likelihood or EM algorithm is a commonly used fitting algorithm for PH distributions. The algorithm is based on the hyper-Erlang distribution (HErD), as shown in Fig. 2. As shown in Fig. 2, the HErD contains M Erlang branches. The initial probability, number of phases (states), and scale parameter of the m-th Erlang distribution are respectively α m (1 ≤ m ≤ M , α m ≥ 0, α 1 + α 2 + · · · + α m = 1), r m , and λ m . The pdf and cdf of the HErD are as follows: Once α = [α 1 , α 2 , · · · , α M ] , λ = [λ 1 , λ 2 , · · · , λ M ] and r = [r 1 , r 2 , · · · , r M ] are determined, the PH distribution is completely defined. For a given M and N , the EM algorithm determining α and λ is described in Table 2 below [24].
Below is a description for determining r. Given the total number of states (i.e., the order) N = r 1 + r 2 + · · · + r M FIGURE 2. State transition graph of an HErD. and the number of possible branches M, all possible situations R = R 1 , R 2 , · · ·, R S N can be obtained, with S N being the number of possible cases of state transition for the set (N, M) and R i = [r 1 , r 2 , · · ·,r M ] being the number of states of each branch in the i-th case. The integer splitting method is used to obtain S N . For example, if N = 5 and M = 1, 2, 3, then S N = 5; i.e., there are five cases of [5], [1,4], [2,3], [1,1,3] and [2,2,1].
The values of N and M are determined based on a trial method. Set the appropriate N and M , and find the best parameter Θ * according to Table 2. Thereafter, the results of PH fitting using the optimal parameters are compared with empirical results. If the fitting is not satisfactory, N and M should be increased. For larger N and M , however, the number of possible cases of the state transition S N grows exponentially. Thus, for large values of N , i.e., N > 10, it is recommended that the HErD is fit with at most three Erlang branches [17]. A larger N would produce better results, but the amount of calculation might be larger, so in the example of this paper, N = 15 and M = 1, 2, 3 is considered sufficient to satisfy the fitting accuracy.

D. EXAMPLE OF PH FITTING
An example of fitting a lognormal distribution is given below. Assuming that a lognormal distribution L has µ = 2 and σ = 0.5, first generate K samples of lognormal distribution L with MATLAB software, and then use the method in Part C of this section to draw the pdf graph. Finally, this graph is compared with the lognormal distribution pdf graph generated by MATLAB.
When using the PH distribution fit, set N = 15 and M = 1, 2, 3. The comparison results are shown in Fig. 3 below.
The result of PH fitting is close enough to the true pdf curve when the sample size is 500; the fitting curve and the real curve are essentially coincident when the sample size is 2000.  to a certain extent [29] and the Global Hydrology Resource Center website [30] can be used to easily obtain the optical radiation energy values of different time zones in different regions, this parameter is selected to characterize the intensity of the lightning in this paper.

III. SHOCK AND DAMAGE MODELS UNDER LIGHTNING
We downloaded lightning data from Chongqing in 2018 from the Global Hydrology Resource Center. As an example, the data from 2018 are used to represent the overall lightning situation in Chongqing, and the optical radiation energy value is selected to characterize the shock magnitude [30]. The specific data are shown in Table 3.
There are 491 data points in Table 3. According to the fitting case in Section II, it is sufficient to approximate the true distribution fitting with 500 data points. Therefore, the algorithm of Section II can be used to fit the opti-  (7), as shown at the bottom of the next page.
The fitting results are shown in Fig. 4. The fitting results are somewhat similar to those of the lognormal distribution. In [31], following analysis of the data of multiple regions, it was found that the light radiation of lightning approximately follows the lognormal distribution. It is evident from this aspect that the results obtained by the above PH fitting method are credible.

2) MODELING OF ARRIVAL TIMES
The shock arrival time refers to the time when lightning occurs and strikes the smart meter. Because the occurrence time of lightning is hard to express with a common distribution and because the probability of striking a smart meter is difficult to describe, it may be assumed that the arrival time obeys a Poisson distribution with a parameter of λ * . In addition, to better represent reality, we divide shocks into three types according to their magnitude.
VOLUME 8, 2020 Suppose there are two shock thresholds D 0 and D 1 ; then, the random shocks whose magnitudes are in the intervals [0, D 0 ], [D 0 , D 1 ], and [D 1 , ∞] are classified into the safety zone, the damage zone and the fatal zone, respectively. A shock in the damage zone generates cumulative damage to the system; a shock in the fatal zone causes immediate failure of the system, i.e., hard failure; and a shock in the safety zone has no effect on the system's failure behavior. According to [32], if the homogeneous Poisson process with intensity λ * is randomly divided into several subprocesses with probability P 1 , P 2 . . . P n , where P 1 + P 2 + . . . + P n = 1, then the obtained subprocesses are independent homogeneous Poisson processes with intensities P 1 λ * , P 2 λ * . . . P n λ * . Let P 1 , P 2 , and P 3 denote the probability that the i-th shock belongs to the damage zone, the fatal zone and the safety zone, respectively, where P 1 + P 2 + P 3 = 1. Assuming that the magnitude of each shock is an i.i.d. random variable, the values of P 1 , P 2 and P 3 can be easily determined from W i . Then, the intensities of the three types of shock can be obtained.
It is well known that if the number of occurrences of a shock obeys the Poisson distribution, the time interval at which the shock occurs is subject to an exponential distribution. Let us define N 1 (t) , N 2 (t) , and N 3 (t) as the number of shocks that have fallen into the damage, fatal and safety zones prior to time t, respectively. If Z k is used to represent the time between two shocks, then Z k ∼ PH (η k , Z k ) .T m = m k=0 Z k ∼ PH (t m , T m ), where t m and T m can be expressed as follows: P(N (t) = m) can be expressed by the following formula: Assume D 0 = 67, D 1 = 1577, and λ * = 0.05. According to the pdf of the shock magnitude and shock thresholds D 1 and D 2 , the lightning arrival time of the three zones obeys the Poisson distribution with intensities of 0.005, 0.045, and 0.005. Taking the shock in the damage zone as an example, the shock arrival time Z k ∼E(0.045), and the PH parameter of the shock arrival time can be expressed as follows: When the m shock arrives, its PH parameter can be expressed as follows: Finally, using formula (8), the probability P(N 1 (t) = m) can be obtained for m shocks in the damage zone before time t. For the same reason, it is also possible to know the probability of m shocks in the safety zone and the fatal zone before time t.

B. ANALYSIS OF SHOCK DAMAGE
When lightning in a damaged area strikes the smart meter, it will cause damage to the smart meter. If the smart meter does not fail after these strikes, the system should satisfy m j=0 Y j < H (12) where Y j represents the damage caused by the j-th shock and H denotes the failure threshold of performance degradation. It is assumed that Y j ∼PH (γ j , Y j ) with order n Y j ; then, S m = m j=0 Y j ∼PH (s m , S m ) with order n Y 1 + n Y 2 + · · · + n Y m . The expressions of s m and S m are given by We assume that lightning acting on a smart meter in the damage area is equivalent to overvoltage, and we carry out overvoltage stress shock tests (data from [33]) to study the damage changes of smart meters under overvoltage. The basic error results for each overvoltage stress are shown in Table 4.
The basic error measured for each sample in Table 4 under overvoltage is shown in Fig. 5.
Because the amount of data per group is small, which may cause the error fitting with PH to be large, we consider whether the data of Table 4 obey a normal distribution. This test is performed using the AD statistic in the software Minitab. The results of the goodness of fit test are shown in Fig. 6.
In the figure, P is the value obtained by the AD detection method. When the P value is greater than 0.05, it is considered to support the hypothesis. The P value in Fig. 6 shows that the data obey a normal distribution under all stress levels when the significance level is 0.05.
To facilitate subsequent calculations, we should consider the relationship among the shock damage under different stresses. Therefore, we need to verify whether the data as a whole are subject to a normal distribution. Hypothesis tests of the mean difference and the variance ratio for the normal population under each stress are conducted with the t distribution and the F distribution, respectively. Table 5 shows the hypothesis test results of the variance ratio, and Table 6 shows the hypothesis test results of the mean difference.
When the significance level is 0.05, the receiving domain is [0.248, 4.03]. The test results in Table 5 are compared with the   receiving domain, and it can be considered that the variance at the eight stress levels is the same when the significance level is 0.05. When the significance level is 0.01, the rejection VOLUME 8, 2020  domain is |T | > 2.8784. The test results in Table 6 are compared with the rejection domain. It can be considered that the mean values at the eight stress levels are the same when the significance level is 0.01.
After the above calculation, it can be considered that under a certain overvoltage stress, the basic error of the smart meter has nothing to do with the magnitude of the overstress and approximately obeys a normal distribution. The mean and standard deviation of the normal distribution can be further determined as Thus, Y j ∼N (0.03, 0.02). Because the result of the PH fitting is a positive distribution and there are some negative parts in Y j , we need to add a value A to all the data to ensure that the value falling in the negative part is small and then sample for PH fitting. Therefore, m × A needs to be considered when referring to the soft failure threshold, where m is the number of times a shock arrives in the damage zone.  (14), as shown at the bottom of the next page.
The fitting results are shown in Fig. 7. (The pdf graph of A + Y j generated by MATLAB is used as a comparison.)

C. SHOCK MODEL UNDER LIGHTNING
After the above analysis, the reliability evaluation of smart meters considering lightning strikes can be expressed as where N 1 (t) and N 2 (t) are defined as the number of shocks in the damage zone and the fatal zone, respectively; Y j represents the damage caused by the jth shock; and H indicates the failure threshold of performance degradation (soft failure threshold).
When obtaining the parameters of the PH distribution of W i ,Y j , N 1 (t) and N 2 (t), the reliability evaluation of the smart meter considering the lightning shock is performed by using (15). Table 7 gives the relevant parameters. Fig. 8 shows

IV. DEGRADATION-SHOCK MODEL UNDER LIGHTNING A. RELIABILITY MODELING IN THE DEGRADATION-SHOCK PROCESS
Considering that the performance of a smart meter degrades under the influence of environmental stress, in addition to the assumptions in Section III, new assumptions must be included in the degradation-shock model. The hypotheses of the degradation-shock model under lightning shock are as follows: 1) The degradation-shock model includes soft damage and hard failure. When the soft damage surpasses the soft failure threshold (as shown by line H in Fig. 9(a)) or the shock magnitude reaches the hard failure threshold (as shown by D 1 in Fig. 9(b)), the smart meter fails.
2) The performance degradation of smart meters consists of two parts: One is the amount of degradation caused by environmental stress, and the other is the amount of degradation that is suddenly increased by each external shock. The total degradation of smart meters is the sum of the two parts.
3) When the smart meter is subjected to a high-intensity lightning strike, as shown by W 4 in Fig. 9(b), it fails; when the smart meter is hit by a medium-intensity lightning strike, such as W 1 and W 3 in Fig. 9(b), the damage can increase, such as Y 1 and Y 2 in Fig. 9(a); when the smart meter suffers a low-intensity lightning strike, the meter is not affected.
Therefore, the reliability of the smart meter's normal operation can be expressed as i=0 X i represents the performance degradation of the smart meter caused by environmental stress, and X i is the amount of performance degradation within t(usually, t = 1).
It is assumed that VOLUME 8, 2020 n X 2 + · · · +n X t/ t . x t and X t can be expressed as Because X i and Y j are subject to the PH distribution, D m also obeys this distribution. If D m ∼ PH (d m , D m ) with order n D m = n X t + n S m , the expressions of d m and D m are given by where X t · 1 + x t = 0. According to (2), the following expression holds: With (9), (16) and (19), the reliability of smart meters under degradation-shock loads can be evaluated using PH fitting.

B. ANALYSIS OF THE DEGRADATION PROCESS CAUSED BY ENVIRONMENTAL STRESS
We consider a Wiener process based on an accelerated stress test to determine X i . This paper selects the linear drift Wiener process to characterize the degradation process: where X (t) is the degradation value of time t; B(t) ∼ N (0, t); µ is the drift coefficient, which reflects the degradation rate; and σ is the diffusion parameter, which indicates the individual difference of the products under different operating conditions. With (19), the damage increment X t = X (t + 1) − X (t) can be calculated, i.e., X i .
The data of the accelerated degradation test and the modeling method can be found in [33]. The temperature and humidity are selected as 20 • C and 60%, respectively, corresponding to the normal state of smart meter operation. Then, µ and σ can be calculated as µ = 5.35e-06 and σ = 0.003, respectively. To reduce the calculation in the degradation-shock model, t is taken as one year; that is, in the subsequent analysis of this paper, X t represents the damage when the smart meter is subjected to temperature and humidity stresses   for one year. Then, in combination with (21), the result of X t can be calculated, i.e., X i .   PH distribution is as follows (22), as shown at the bottom of the this page.
The fitting results are shown in Fig. 10 (the pdf graph of B + X i generated by MATLAB is used as a comparison):

C. RESULTS AND DISCUSSION
Based on the analysis in Section III and combined with the degradation analysis caused by the environmental stress in Part B of this section, the reliability evaluation results of a smart meter under degradation-shock loads can be obtained by using (16). To verify the accuracy of the results, the Monte Carlo method is used for comparison. The reliability assessment results considering the environment-lightning effect are shown in Figs. 11-14. Fig. 11 shows that the reliability evaluation method proposed in this paper and the Monte Carlo method are essentially consistent, which verifies the accuracy of the reliability evaluation method in this paper. In addition, the reliability curve considering lightning strikes is significantly different from that without considering lightning strikes. Therefore, it is necessary to consider lightning strikes when modeling. When the shock arrival temporal spacing or the shock hard failure threshold is reduced, the reliability of smart meters decreases rapidly, so it is necessary to take effective lightning protection measures in places where lightning is frequent.

V. CONCLUSION
In this paper, a reliability evaluation method for smart meters based on the degradation-shock of a PH distribution is proposed. The method takes into account both environmental stresses (temperature and humidity) and lightning strikes. For some unconventional distributed data processing, the typical approach is to make approximating assumptions, which may reduce the applicability of the model. The PH distribution has favorable characteristics and can approximate any positive continuous distribution, so this paper chooses a model based on this distribution. Finally, through comparison with the Monte Carlo method, the accuracy of the reliability evaluation in this paper is verified.
The reliability calculated based on the method of this paper can provide support for a company in arranging maintenance and replacement decisions to make the best use of smart meters and avoid waste of manpower. Based on the presented results and discussion, the following main contributions can be drawn from our investigation.
(1) This paper proposes a degradation-shock modeling and reliability assessment method based on the PH distribution, which can easily obtain the reliability assessment of smart meters in different regions; that is, the proposed method is universal.
(2) In the reliability assessment model of smart meters, the impact of hard and soft damage caused by lightning and overvoltage on the reliability of the meter is considered. The shock damage and the shock magnitude are nonlinear.
(3) A correlation model of shock soft damage and shock magnitude is established through a smart meter overvoltage stress shock test.
Furthermore, to more comprehensively consider the shock process, additional work can be done in the future. First, the shock arrival time as modeled by the nonhomogeneous Poisson process will be further considered. Second, the safety threshold and damage threshold of the shock process may change with the degree of shock and degradation damage, which will also be considered. In addition, there may be a dependent relationship between the shock damage and degradation damage: The magnitude of the shock damage is related to the current degradation degree.