Mapping Reliability Predictors of Low-Voltage Metal Oxide Surge Arresters Using Contour Plots

Reliability assessment of MOSA devices requires accurate estimation of life distribution parameters or reliability predictors. Estimated reliability predictors consist of reasonable indicators of failure probability and mean life trends as a result of applied stress over a given period of time. Recent literature suggests finite-value or point-based estimation approach of reliability predictors in the context of continuously-applied distorted voltage stress to MOSA devices. However, this technique is prone to inaccuracies, and may not effectively reflect the impact of applied stress on reliability predictors. In this work, the use of contour plots is proposed as an alternative method for discrete mapping of estimated reliability predictors. Therefore, accelerated life tests – consisting of constant electro-thermal stress of: 85% of the reference voltage with and without harmonics, at the temperature of $135 ^{o}C$ for the time-period of 96 hours – are conducted on a set of low-voltage MOSA sourced from different manufacturers. The experimentally obtained times to failure are consistent with the two-parameter Weibull life distribution. The log-likelihood estimates on the Weibull probability density function are invoked using computational algorithms which yield the contour plots. These plots depict the reliability predictor mappings in the $\alpha \beta $ plane. Results show that harmonic content in the applied stress causes extended peak amplitude which represents the point of highest probability of reduced reliability, and thus the highest accuracy point for the estimated shape parameter.


B l
probability density function at (t F ) l (t F ) i measured time to failure t t test time r number of failed samples V 1 mAac ac reference voltage V 1 mAdc dc reference voltage V 1 mAdc change in dc reference voltagê (t F ) i extrapolated time to failure T t test temperature F (i, n) cumulative failure probability f (t F ; β, α) Weibull probability density function i failure ranking The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney. n total number of samples γ (x i , y i ) correlation factor x i logarithmic expression of F (i, n) y i logarithmic expression of(t F ) i x average value of x i y average value of y i β shape parameter β estimated shape parameter β cont extracted shape value within contour β j estimated shape parameter at point j α scale parameter α estimated scale function α cont extracted scale value within the contour α k estimated scale parameter at point k L estimated log-likelihood L cont estimated log-likelihood within contour VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ L jk estimated log-likelihood at α k and β j L (β, α) log-likelihood function r i=1 product of the probability density function χ 2 0.9;1 chi square function at 90% confidence level

I. INTRODUCTION
The continued ac or dc field-induced conduction is reported to be one of the major causes of electrical degradation or reduced reliability of metal oxide varistor (MOV)-based surge arresters [1]- [4]. This long-term failure or progressive reliability reduction process of varistor arresters is commonly experienced in practice [5], [6], and may expose the system to the risk of critical damage, prolonged stalling in production and considerable losses in revenue. To ensure reliable surge protection in power systems, electronic and data circuits: condition monitoring and reliability assessment of MOSA devices is essential practice to be adhered to [7], [8]. Besides thermal stress having proven to be an important non-electrical contributing factor to the reliability reduction, power system harmonics have been reported to be significant aggravating factor to reliability reduction phenomenon of metal oxide surge arresters (MOSA) [9]- [11]. A mathematical model proposed in [12], predicts life acceleration, and hence reliability reduction of MOSAs for a given percentage content of harmonic-distortion in the voltage stress. Therefore, in the context of reliability assessment of MOSA devices, the knowledge of the most suitable life distribution − which best describes the behaviour or performance of MOSA devices under defined operating conditions (electrical and non-electrical stresses) − consists of an essential adopted step. Subsequent to this step, the estimation of reliability predictors of the life distribution is logical, and thus critical to the accuracy of the assessment of reliability. However, estimation techniques of reliability predictors such as currently proposed in the literature mostly suggest a finite value-based approach to determine these predictors [13]. This approach is prone to innacuracies and incorrect reliability assessment. Another approach suggests the determination of the maximum value among a set of possible estimated values. The drawback in the latter approach resides in the possibility of common predictor values for a change in the stress condition applied to similar samples, and therefore may yield to erroneous inferences of reliability responses. This work is aimed at mapping discretely estimated reliability predictors using contour plots. This approach improves the existing estimation techniques and eliminates the possibility of inaccuracies associated with finite value-based estimation of predictors. Furthermore, the responses of discretely estimated reliability predictors to variables or additional factors in the applied stress are described in terms of shape mutation of the obtained maps.
Therefore, three groups of low voltage MOSAs, with similar electrical and physical properties − sourced from three different manufacturers − are subjected to continuously-applied ac field and thermal stresses for a constant time-period. A two-fold ac field stress was applied -i.e. with and without harmonic-distortion. The data (failure times) produced in each of the six test cases conducted fit the two-parameter Weibull statistical distribution. The maximum likelihood estimation (MLE) was applied in order to determine the distribution parameters. The contour plots obtained provide comparative maps of reliability predictors for sets of MOSA devices subjected to different intensities of harmonics in the ac voltage stress. The results indicate at 90 % confidence level that: voltage harmonics cause the peak amplitude of the predictor maps to extend upward and its area to shrink and/or shift leftside. Both the peak amplitude extension and the shifting or shrinking of the predictor maps indicate reduced reliability of MOSA samples when the electric field is distorted with harmonics.

II. A REVIEW OF RELIABILITY PREDICTION TECHNIQUES ON MOSAs
Over the last few decades, reliability prediction studies such as conducted on metal oxide surge arresters and other solid insulating materials are associated with the following fundamental steps: degradation modelling (accelerated life tests) and reliability assessment (life distribution and estimation of reliability predictors): Herr et al. [14] used both types of accelerated life tests (progressive or step-stress and constant stress) in order to generate reliability prediction or assessment models for a signal diode, signal and power transistors, silicon controlled rectifier (SCR), and metal oxide varistor (MOV). The constant stress was found to be suitable for modelling the time relationship between the applied stress and the device performance. The times to failure that constituted the data obtained could be used to build reliability prediction models for both varistor under ac or dc bias combined with temperature which is a non-electrical stress. Cygan and Laghari [15] studied the techniques and models commonly applied in the analysis and prediction of the lifetime of solid insulators (including arresters) when subjected to constant electrical and thermal stresses. This research observed that life tests at higher than normal bias voltage and temperature form the most popular technique to simulate electrical ageing phenomenon. Furthermore, the Weibull and log-normal statistical distributions are reported to be effective in the analysis of life deterioration of solid insulating materials. Escobar et al. [16] described a relationship between degradation strength and component failure. This study concludes that accelerated degradation enables inferences to be made on data and the resulting lifetime distributions. It further suggests that an acceleration model is required to relate the effect of an accelerating variable (e.g. temperature, etc. . . ) to the rate of failure of a component or product, and the maximum likelihood estimation (MLE) is usually relied upon in order to estimate the parameters of the model. Meeker and Escobar [17] formulated a guide to degradation analysis for several types of life tests. This guide suggested that degradation analysis could potentially be used for reliability improvement of products or components. Bae et al. [18] studied the relationship between degradation models (additive and multiplicative with single random effects) and the resulting lifetime distributions. This investigation suggests that despite the direct link between the degradation model and the lifetime distribution, careful attention should be devoted to random effects on the degradation model. The IEEE 930 TM [19] proposed a number of guidelines relative to the application of the two or three-parameter Weibull statistical distribution in the modelling of insulation life data (times to failure). Among the proposed guidelines, the following could be retained: the percentage cumulative degradation of the Weibull cumulative density function (CDF) is obtained using the Ross approximation, the adequacy or the goodness of fit of the CDF obtained is evaluated using the correlation factor. For large data the shape and scale parameters of the Weibull distribution are estimated using the least squares regression method. Meshkatoddini [20] applied the Monte Carlo algorithm for the purposes of statistically analysing the conduction or ageing phenomenon in thin ZnO varistors. He concluded that the number of conducting ZnO grains that provide the current path between the electrodes of thin varistors fits the log-normal statistical distribution. Such an increase in the number of ZnO conducting paths of varistor may be associated to the probability of failure of MOSA devices. Brown [21] studied the ageing phenomena of MOSA devices and observed that the lifetime of an MOV-based arresters could be modelled using the Arrhenius rate equation. The failure rate of MOSAs under normal operation can be estimated using the Arrhenius model if the most significant stress is thermal and the expected mean life (ML) is logarithmically related to the inverse of temperature. Chu and Ke [22] examined the two commonly used computational techniques for the estimation of the Weibull parameters: the MLE and the least squares method (LSM). This study recommends the use of the least squares method when the sample size is small, and that of the MLE when the sample size is large. This observation is supported by the work presented in [23], which investigated the accuracy of these two computational techniques such as applied in the life modelling of MOSAs. Yang et al. [24] employed step-stress or progressive accelerated degradation test as well as multivariate degradation modeling in order to estimate the reliability of rotary encoders. This study supports the pertinency of accelerated degradation tests prior to failure modeling such as applied to MOSAs in [12]. Sun et al. [25] proposed a computer-aided approach, which relies on the physics of failure, in a bid to perform reliability analysis of electronic components. This technique advocates, among other things: the accurate interpretation of multi-correlation that may relate product life (time to failure) or performance reliability with several key parameters or indicators. Therefore, reliability prediction or modelling of failure data is essentially reliant on accurate estimation of the predictors associated to the lifetime distribution resulting from constantly or progressively applied accelerated degradation or life tests to which the samples are subjected to.

III. METHODOLOGY A. ACCELERATED LIFE TESTS
Standard accelerated life tests are usually employed to simulate electro-thermal (eletrical and non-electrical) degradation of MOSAs and/or other solid insulators [26], [27]. For the purposes of this test, the Nabertherm P 330 heat chamber was configured to provide constant environmental temperature of 135 o C. The continuous basic wave of the applied voltage stress was obtained from the TDGC 2 50 Hz ac variable voltage source, and amounted to 85% of the reference voltage (V 1 mAac ) of the 360 MOSA samples used. Each manufacturer group consisted of two sub-groups of 60 samples. As an acceleration degradation factor, odd harmonic-frequencies 3 rd , 5 th , 7 th were induced and superimposed on the basic wave using a triac-based ac voltage controller switching a resistive-inductive load. The ac voltage controller unit as well its load was connected across the MOSA under test. For each test run, 5 samples were connected in parallel or across each between two terminals inside the heat chamber and supplied from the voltage source. The life test time consisted of 96 hours of continously applied electro-thermal stresses. The connection between the ac voltage supply system to the terminals consisted of single-core, 1.5 mm 2 Silflex high temperature conductor. For the purposes of time-synchronisation between electrical and heat stresses, the voltage source was time-delayed until thermal transients had subsided to fixed test temperature. The two-channel MT250 and the three-channel K5020 data loggers, which essentially record and provide patterns of voltage across MOSAs each 30 seconds of time, were employed for capturing times to failure data. The experimental set-up is indicated in Figure 1.

B. DESCRIPTION OF ARRESTER SAMPLES
The tested MOSA samples consisted of low voltage units of 20 mm size, and sourced from three manufacturers: A, B and C. Therefore, A 1 , B 1 and C 1 will refer to samples from manufacturer A, B and C, tested without harmonics in the voltage stress. Similarly, A 2 , B 2 and C 2 will represent samples from respective manufacturers and tested with harmonics in the voltage stress. In order to ensure accurate analysis of the effect of voltage harmonics on the reliability predictors, life distribution resulting from A 2 samples should be compared to that of A 1 , B 2 with B 1 and C 2 with C 1 . The physical properties of the samples used, and the magnitudes of applied electrical stresses, are provided in Table 1.

C. FAILURE CRITERION, DATA FITNESS AND ESTIMATION OF WEIBULL PARAMETERS 1) FAILURE CRITERION
For the purposes of this work, an arrester sample is deemed to have failed provided that the recorded time to failure is less than the test time and the percentage difference between the V 1 mAdc points of the V − I curve, measured before and after accelerated life test, is higher or equal to 5%. This could be expressed as follows: The times to failure data obtained from accelerated life tests are extrapolated using the following equation [23], [26]: For each test conducted, the number of failed MOSA samples is indicated in Table 2. Prior to statistical inferences, the fitness of the obtained data (times to failure) to the Weibull distribution should be tested [28], [29]. This requires the correlation factor to be determined. Therefore, the Ross approximation is applied in order to determine the cumulative failure probability.
The correlation factor is obtained using the following relationships: The logarithmic expressions of the cumulative failure probability and that of the time to failure yield the following equations: and Similarly, the average values of the logarithmic expressions are given as follows: x = x i r (6) and y = y i r (7) Parametric analysis is conducted as the actual failure time quantities are employed and the underlying distribution is known. Furthermore, for data fitness, the correlation factor must be higher or equal to the critical coefficient value provided in [28], [30].

3) ESTIMATION OF THE WEIBULL PARAMETERS
For the purpose of the data available in this work, the application of the LSM in order to estimate the shape and scale parameters of the Weibull is followed. The LSM is a widely used approach in estimating Weibull parameters. The validity of the point estimation of predictors or parameters obtained in this process can subsequently be tested using contour plots. Therefore, the shape and scale parameters can be determined using equations (8) and (9) such as proposed in [29]: Equations (8) and (9) yield finite values of the shape and scale parameter, respectively. Such values only represent estimates within a wider range of possiblities. The correlation factors as well as the finite values of the estimated parameters using LSM are provided in Table 3. Since the reliability function is dependant on both the shape and scale parameters, reliance on a wider range or interval of estimated parameters may provide a clearer indication of the upper and lower bounds within which the predictors could be defined. Ultimately, the behaviour of the reliability function could be observed over a range rather than on the basis point estimates. Instead, the likelihood ratio bounds is applied here to determine the extreme values of the contour plots on both axes for a given confidence level [31]- [33].

IV. ALGORITHMS FOR MAPPING DOMAINS OF RELIABILITY PREDICTORS
This set of algorithms enables mappings of the failure density and, ultimately, of the reliability function in the αβ plane to be produced. The Weibull distribution of the times to failure data, such as discussed in section III, prompts the log-likelihood function [L(β, α)] to be invoked in terms of the Weibull probability density function [f (t F;α,β) )]. The log-likelihood function is expressed in terms of the following equation: where: The log-likelihood estimates − corresponding to the scale and shape parameters (predictors) of the Weibull life distribution − are computed using (10). In order to determine the contour plots, Matlab software was used to implement the presented algorithms for mapping of the domains of the reliability predictors and to perform the analyses. The following set of algorithms are developed:

Algorithm 1 Estimation of the Log-Likelihood Matrix
Input t F Initialize α, β ranges α k , β j for j = 1 to length of β j do for k = 1 to length of α k do for l = 1 to length of t F do

end for end for
Through algorithms 1 and 2,L and L jk are obtained.L corresponds to multiple α and β, other thanα andβ, that naturally forms a closed contour shape on the αβ plane. Since these are discrete estimates, a 'band' or interval about L is required to obtain a suitable contour. For this reason, L jk is computed as it contains the log-likelihood estimates corresponding to the selected ranges of α and β. The 'band' or Algorithm 2 Extraction ofL and L jk to Obtain Contour Input α k , β j , L jk for j = 1 to length of β j do for k = 1 to length of α k do if L jk ∈ [0.7L; 1.3L] then L cont ← L jk α cont ← α k β cont ← β j end if end for end for interval, aboutL, can now be extracted from L jk . This is done using algorithm 2, which extracts values of L jk that lie within the pre-determined interval -i.e. L cont . Thereafter, the α k and β j that yield L cont values within the range are extracted -i.e. α cont and β cont -and used to plot the contour. For this work, an interval of [0.7L; 1.3L] was selected.

V. RESULTS AND ANALYSIS
The resulting contour plots or predictor mappings obtained are indicated in figures 2 (A 1 and A 2 samples), 3 (B 1 and B 2 samples) and 4 (C 1 and C 2 samples). Detailed examination of these figures reveal the following observations: 1) The peak amplitude of the predictor mappings, obtained from samples tested with distorted ac field stress, appears to be relatively higher than that of samples treated without distortion. For the three cases observed, this suggests that the maximum value of the shape parameter is located in contour plots A 2 , B 2 and C 2 , respectively. This could be summarized as follows:β max for A 1 samples <β max for A 2 β max for B 1 samples <β max for B 2 β max for C 1 samples <β max for C 2 2) The base width of the predictor mappings, obtained from samples tested with distorted ac field stress, prove relatively shorter than that of samples treated without distortion in the ac field stress. This implies that the maximum value of the scale parameter is located in contour plots: A 1 , B 1 and C 1 , respectively. This could be expressed as follows: α max for A 1 samples >α max for A 2 α max for B 1 samples >α max for B 2 α max for C 1 samples >α max for C 2 3) The common area between the predictor mappings obtained in each of the cases observed indicates a set of equal estimated predictors to both life distributions despite differences in terms of applied stress. This consists of an important observation which highlights how the contour plot-based method improves the previously used techniques of predictor estimation. VOLUME 8, 2020  4) The contour plots obtained in the three cases observed have proven to be distinct in shape. This is owed to manufacturers' differences in microstructure additives of the MOSA samples tested. The change in shape of the contour plots obtained is indicative of the response of the predictors to THD content in the applied stress. 5) It could also be observed that the point-based estimated parameters such as given in Table 3 are contained in the contour plots obtained, and are in some cases located in the common area. 6) The peak amplitude extension as well as the width shrinkage or shift of the contour plots, such as observed in all cases, are consistent with the THD content in the applied voltage stress. Since the peak amplitude consists the highest probability of reduced reliability, it could be regarded as the accuracy point as far as one of the predictors (the shape) is concerned.
Since the contour plots, such as obtained in this work, result from the log-likelihood of the Weibull probability density function, whereby the peak amplitude of the curve represents the highest value of the shape parameter (probability of reduced reliability), and the length between the two sides of the curve is indicative of the highest value of the scale parameter (mean life) of the distribution. Therefore, both the upward extension in terms of height and the leftside shift and/or the shrinkage of the area under the contour plot consist of key visual indicators, used in this approach, in order to evaluate the reliability of a distribution. Therefore, these observations infer higher probability of failure and shorter mean life for MOSA distributions exposed to higher THD content of harmonics in the applied stress. This implies that A 2 , B 2 and C 2 samples have proven to be less reliable than A 1 , B 1 and C 1 , respectively.

VI. CONCLUSION
Contour plots approach provide an alternative technique to effectively determine the mappings of discretely estimated predictors. This is essential for reliability assessment or evaluation of life distributions for MOSAs or other solid insulators. The closed contour shape on the αβ plane, such as obtained through computational algorithm, provides multiple estimates of the predictors. This method overcomes the possibility of inconclusive outcome that may result from comparative evaluation of reliability or life distributions with equal finite value or point-based estimated predictors. The shape and area of the contour plots constitute key indicators of predictor response to harmonics in the applied stress. Furthermore, this technique enables the determination on how close or accurate predictors obtained from point-value estimation technique may be to the peak amplitude of the contour plot.