# A Poisson-Fault Model for Testing Power Transformers in Service

Authors: Dengfu Zhao; Zheng Zhao; Qihong Duan; Gongnan Xie
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper presents a method for assessing the instant failure rate of a power transformer under different working conditions. The method can be applied to a dataset of a power transformer under periodic inspections and maintenance. We use a Poisson-fault model to describe failures of a power transformer. When investigating a Bayes estimate of the instant failure rate under the model, we find that complexities of a classical method and a Monte Carlo simulation are unacceptable. Through establishing a new filtered estimate of Poisson process observations, we propose a quick algorithm of the Bayes estimate of the instant failure rate. The proposed algorithm is tested by simulation datasets of a power transformer. For these datasets, the proposed estimators of parameters of the model have better performance than other estimators. The simulation results reveal the suggested algorithms are quickest among three candidates.


---

## Body

## 1. Introduction

The failure of power transformers is one of most frequent causes of long interruptions of power systems. Many researchers have studied evaluation schemes of a piece of equipment of power systems; see [1–3] and reference therein, and the IEEE C57-110-98 standard [4] is developed to determine the capability of a power transformer. However, due to the long lifetime of a power transformer as well as the complexity of an experiment, few statistical experiments about the lifetime or the failure rate of power transformers have been performed. A singular experiment, in which eight large transformers are tested, is reported in [5]. On the other hand, there are many power transformers in service, and we can obtain data of periodic inspections and maintenance for them. The fact inspires us to find from the data the relation between the failure rate and the working condition of a power transformer.While investigating a statistical method for this problem, we face three difficulties. The first difficulty is as follows: data is collected according to a standard test procedure [6, 7], but the procedure does not take into account our problem. What we know is the working condition of a power transformer at an inspection as well as whether or not a transformer power fails. As time between successive inspections is far less than a typical lifetime of a transformer power, we can assume that the number of failures that occur between two successive inspections has a Poisson distribution; see [8, 9] for more details about a Poisson model of similar cases. The rate of the Poisson distribution, which is proportional to the instant failure rate, is associated with the current working condition of a power transformer. According to an idea of Bayesian experimental design [10], we assume that, for a power transformer, the instant failure rate with a given working condition has a prior normal distribution. Then the Bayes estimate is the best point estimate to the instant failure rate with the quadratic loss function.The second difficulty is that the mean and the variance matrix of the prior distribution might not be easily determined. Some empirical models, including the inverse power law model [11] and the Arrhenius model recommended by the IEC loading guide [12], may motivate us. These models study the influence of working temperature and operating environment on the lifetime of a power transformer. Substituting in each model by an appropriate parameter, we find that the model is an ideal case of an autoregression model. Considering observed errors as well as possible bedrock errors of these empirical models, we assume that the instant failure rate of a power transformer with respect to working conditions is an autoregressive process. Furthermore, by applying a classical method of time series and a technique of martingales, we can establish estimates to parameters for the prior distribution. Comparing with [6], we establish an acceptable method for investigating the instant failure rate of a power transformer.The last difficulty appears while computing the Bayes estimate to the instant failure at different working conditions. LetN denote the total number of working conditions and let L denote the number of failures. The Bayes estimate involves a Lth moment of N-dimensional normal random variable. The time complexity of a classical algorithm to the moment is O(NL). Although a Monte Carlo simulation is not sensitive to the dimension N, its performance is not good when L is slightly large.Recently, there is some advance in filtered estimates for Poisson process observations; see [13, 14] and references therein for example. By modulating the preceding estimates for our case, we can establish an evolutionary equation for the Bayes estimate. A closed-form solution to the evolutionary equation is obtained and a corresponding quick algorithm is then established. Similarly to [8, 13], we establish a closed-form solution to the problem. Moreover, the corresponding algorithm is exact. As a comparison, [13] used an approximation algorithm. Finally, the simulation result reveals that the algorithm is quicker than a Monte Carlo simulation, which is suggested by [15] for a similar problem.The rest of the paper is organized as follows. We construct the model in Section2 and discuss its relation with available empirical models. Section 3 presents the setting of the prior distribution. In Section 4, we investigate the Bayes estimate of the instant failure rate. After establishing a classical algorithm and a Monte Carlo simulation of the Bayes estimate and exploring its performance, we present a quick algorithm according to a new theory in Section 4. Simulation results are presented in Section 5. Section 6 concludes the paper.


## 2. Constructing a Poisson-Fault Model

In this section, after a literature review about some empirical models, we propose a Poisson-fault model of a power transformer.


### 2.1. Empirical Models

The Arrhenius model is accepted widely. The IEC loading guide [12] uses the model and many authors investigated it; see [6, 7, 16] for more details. According to the Arrhenius reaction rate equation, the average lifetime t of a piece of equipment at different temperature in the Arrhenius model is given by
(1)t=Aexp⁡(-kT).
Here T is the temperature, A is a nonthermal constant, and k is the multiplication of the activation energy and the Boltzmann constant.The power law model [11] can describe dielectric failures of a capacitor under increased voltages and the lifetime of a multicomponent system. The average lifetime t in that model is given by
(2)t=γV-η.
Here V represents a voltage level or a level of a severe operating environment. With the nonlinear transformation U=log⁡(V) of the level, we have that t=γexp⁡(-ηU). Hence, we may regard the power law model as another edition of the Arrhenius model (1).Typically,T in (1) may be substituted by other parameters. In most applications, it is assumed that the increased level of the working condition affects the average of the lifetime and the corresponding instant failure rate.


### 2.2. Poisson-Fault Model

In general, it is assumed that the number of failures has a normal distribution in a given time interval. The assumption is correct if test duration is long. As we know, when using a short test duration, the number of failures has a Poisson distribution; see [8, 9] for more details.We consider a right-continuous and complete filtrationℱ={ℱn:n=0,1,…} on a complete probability space. Two nonnegative stochastic processes X={Xn:n=1,2,…} and M={Mn:n=1,2,…} are adapted to the filtration ℱ. We refer to [17] for an introduction of filtrations and stochastic processes adapted to filtration. We call M a Poisson-fault model with the intensity X if these two processes satisfy the following conditions.(C1)
Forn=1,2,…, Eϵn=0,   Dϵn is a positive constant σ2 and ϵn+1 is independent with respect to ℱn. For given constants a, b such that |a|<1, b>0,
(3)Xn+1=aXn+b+ϵn+1.(C2)
The distribution ofMn+1 conditional on a σ-field ℱn∨σ{Xn+1} is a Poisson distribution with a parameter Xn+1. That is, for k=0,1,…,
(4)P(Mn+1=k∣ℱn∨σ{Xn+1})=Xn+1kk!exp⁡(-Xn+1).Here, ℱn describes the full information at level 1,…,n. And the combination ℱn∨σ{Xn+1} describes the information of both ℱn and σ{Xn+1}. In the above model for a power transformer, n represents the level of working conditions. Mn represents the number of failures between two successive inspections at level n.


### 2.3. Relation between Models

To motivate the model, suppose thatFn is the instant failure rate of a power transformer under a working condition with level n. When there are p power transformers at a working condition and the length of time interval between inspections is d, we have that the average number Xn of failures is proportional to the failure rate Fn. Concretely,
(5)Xn=dpFn.
Moreover, Mn, the number of failures at the working condition n, has a Poisson distribution with parameter Xn. Hence, we conclude that a Poisson-fault model describes an instant failure rate under different working conditions, and an instant failure rate is an autoregressive process.Now we point out the relationship between the Arrhenius model (1) and our model. We consider a special case in which b=0 and ϵn=0 for n=2,3,…. In this case, we have
(6)Xn=an-1X1.
We may write the corresponding “temperature” T=C-nΔ, where C is a constant and Δ is the change of two successive working conditions. Moreover, the failure rate Fn is the reciprocal of the average lifetime for a piece of equipment for a short time interval d. Hence, we have that the average lifetime t at the level n is given by t=Aexp⁡(-kT), where A=a(1-C)/Δ(dpX1)-1 and
(7)k=-log⁡(a)Δ.By comparing this equation with (1), we may conclude that the autoregressive process X incorporates and extends features of the Arrhenius model (1).


## 2.1. Empirical Models

The Arrhenius model is accepted widely. The IEC loading guide [12] uses the model and many authors investigated it; see [6, 7, 16] for more details. According to the Arrhenius reaction rate equation, the average lifetime t of a piece of equipment at different temperature in the Arrhenius model is given by
(1)t=Aexp⁡(-kT).
Here T is the temperature, A is a nonthermal constant, and k is the multiplication of the activation energy and the Boltzmann constant.The power law model [11] can describe dielectric failures of a capacitor under increased voltages and the lifetime of a multicomponent system. The average lifetime t in that model is given by
(2)t=γV-η.
Here V represents a voltage level or a level of a severe operating environment. With the nonlinear transformation U=log⁡(V) of the level, we have that t=γexp⁡(-ηU). Hence, we may regard the power law model as another edition of the Arrhenius model (1).Typically,T in (1) may be substituted by other parameters. In most applications, it is assumed that the increased level of the working condition affects the average of the lifetime and the corresponding instant failure rate.


## 2.2. Poisson-Fault Model

In general, it is assumed that the number of failures has a normal distribution in a given time interval. The assumption is correct if test duration is long. As we know, when using a short test duration, the number of failures has a Poisson distribution; see [8, 9] for more details.We consider a right-continuous and complete filtrationℱ={ℱn:n=0,1,…} on a complete probability space. Two nonnegative stochastic processes X={Xn:n=1,2,…} and M={Mn:n=1,2,…} are adapted to the filtration ℱ. We refer to [17] for an introduction of filtrations and stochastic processes adapted to filtration. We call M a Poisson-fault model with the intensity X if these two processes satisfy the following conditions.(C1)
Forn=1,2,…, Eϵn=0,   Dϵn is a positive constant σ2 and ϵn+1 is independent with respect to ℱn. For given constants a, b such that |a|<1, b>0,
(3)Xn+1=aXn+b+ϵn+1.(C2)
The distribution ofMn+1 conditional on a σ-field ℱn∨σ{Xn+1} is a Poisson distribution with a parameter Xn+1. That is, for k=0,1,…,
(4)P(Mn+1=k∣ℱn∨σ{Xn+1})=Xn+1kk!exp⁡(-Xn+1).Here, ℱn describes the full information at level 1,…,n. And the combination ℱn∨σ{Xn+1} describes the information of both ℱn and σ{Xn+1}. In the above model for a power transformer, n represents the level of working conditions. Mn represents the number of failures between two successive inspections at level n.


## 2.3. Relation between Models

To motivate the model, suppose thatFn is the instant failure rate of a power transformer under a working condition with level n. When there are p power transformers at a working condition and the length of time interval between inspections is d, we have that the average number Xn of failures is proportional to the failure rate Fn. Concretely,
(5)Xn=dpFn.
Moreover, Mn, the number of failures at the working condition n, has a Poisson distribution with parameter Xn. Hence, we conclude that a Poisson-fault model describes an instant failure rate under different working conditions, and an instant failure rate is an autoregressive process.Now we point out the relationship between the Arrhenius model (1) and our model. We consider a special case in which b=0 and ϵn=0 for n=2,3,…. In this case, we have
(6)Xn=an-1X1.
We may write the corresponding “temperature” T=C-nΔ, where C is a constant and Δ is the change of two successive working conditions. Moreover, the failure rate Fn is the reciprocal of the average lifetime for a piece of equipment for a short time interval d. Hence, we have that the average lifetime t at the level n is given by t=Aexp⁡(-kT), where A=a(1-C)/Δ(dpX1)-1 and
(7)k=-log⁡(a)Δ.By comparing this equation with (1), we may conclude that the autoregressive process X incorporates and extends features of the Arrhenius model (1).


## 3. Prior Distribution and Its Parameter

For the latent autoregressive processX, (X1,…,XN) has a prior normal distribution. In this section, we propose a method to estimate parameters of the prior distribution.


### 3.1. Background of Time Series

Here we review some results of the classical theory of time series; see [18] for a detailed introduction to time series.We have that forn≥m, Cov(Xm,Xn+1)=aCov(Xm,Xn) and DXn+1=a2DXn+σ2. Provided that Xn is a stable stochastic process, we have EXn=b/(1-a) and DXn+1=σ2/(1-a2). That is, the parameter of the prior distribution is determined by providing a, b, and σ2. Now we turn to estimate these key parameters.With the conditions|a|<1 and b≥0, the intensity process Xn is an asymptotic stable stochastic process. Hence, three sequences EXn, EXn2, and EXn+1Xn are convergent. Therefore, these sequences are bounded.Assume(x1,…,xN) is a dataset of observations for the intensity process Xn, n=1,…,N. It is well known that we can solve linear equations
(8)A(a^b^)=β,A=(∑n=1N-1‍xn2,∑n=1N-1‍xn∑n=1N-1‍xn,N-1),β=(∑n=1N-1‍xn+1xn∑n=1N-1‍xn+1)
for estimators a^ and b^. As N tends to ∞, a^ and b^ converge to a and b with a probability 1, respectively. Moreover, as Eϵn=0, we can obtain an estimator σ^2 by the following formula:
(9)σ^2=1N-1∑n=1N-1‍(xn+1-a^xn-b^)2.


### 3.2. Parameter Estimation to the Prior Distribution

Assume that we face a history record of failures. We can then obtainMn, n=1,…,N by counting failures. However, we cannot obtain exactly Xn due to noises. In this case, we must establish new estimators for a and b only based on observations of M1,…,MN. In a simple heuristic, we may replace Xi by Mi in (8) to (9) and obtain estimators a* and b* for a and b, respectively. Based on simulation results in Section 5, we will conclude that a* and b* are bad estimators. In the following, we construct a good approximation of Xn based on M1,…,MN and obtain better estimators for a and b.Theorem 1.
M
=
{
M
n
:
n
=
1,2
,
…
} is a Poisson-faults model with an intensity X={Xn:n=1,2,…}. As N tends to ∞, with a probability 1, one has that
(10)1N-1∑n=1N-1‍(Mn-Xn)⟶0,1N-1∑n=1N-1‍(Mn+1Mn-Xn+1Xn)⟶0,1N-1∑n=1N-1‍(Mn2-Mn-Xn2)⟶0.See AppendixA for the proof of Theorem 1.Then, it follows from Theorem1 and (8) that we can solve linear equations
(11)A~(a~b~)=β~,A~=(∑n=1N-1‍Mn2-Mn,∑n=1N-1‍Mn∑n=1N-1‍Mn,N-1),β~=(∑n=1N-1‍Mn+1Mn∑n=1N-1‍Mn+1)
for estimators a~ and b~. Moreover, it follows from Theorem 1 and the convergence of a^ and b^ that a~ and b~ converge to a and b with probability 1, respectively. Therefore, as N tends to ∞, a~ and b~ converge to a and b with a probability 1, respectively. Similarly, from (9), we can construct an estimator
(12)σ~2=b~2+1N-1∑n=1N-1‍a~2(Mn2-Mn)+(Mn+12-Mn+1)-2a~MnMn+1-2b~Mn+1+2a~b~Mn
for the parameter σ2.


## 3.1. Background of Time Series

Here we review some results of the classical theory of time series; see [18] for a detailed introduction to time series.We have that forn≥m, Cov(Xm,Xn+1)=aCov(Xm,Xn) and DXn+1=a2DXn+σ2. Provided that Xn is a stable stochastic process, we have EXn=b/(1-a) and DXn+1=σ2/(1-a2). That is, the parameter of the prior distribution is determined by providing a, b, and σ2. Now we turn to estimate these key parameters.With the conditions|a|<1 and b≥0, the intensity process Xn is an asymptotic stable stochastic process. Hence, three sequences EXn, EXn2, and EXn+1Xn are convergent. Therefore, these sequences are bounded.Assume(x1,…,xN) is a dataset of observations for the intensity process Xn, n=1,…,N. It is well known that we can solve linear equations
(8)A(a^b^)=β,A=(∑n=1N-1‍xn2,∑n=1N-1‍xn∑n=1N-1‍xn,N-1),β=(∑n=1N-1‍xn+1xn∑n=1N-1‍xn+1)
for estimators a^ and b^. As N tends to ∞, a^ and b^ converge to a and b with a probability 1, respectively. Moreover, as Eϵn=0, we can obtain an estimator σ^2 by the following formula:
(9)σ^2=1N-1∑n=1N-1‍(xn+1-a^xn-b^)2.


## 3.2. Parameter Estimation to the Prior Distribution

Assume that we face a history record of failures. We can then obtainMn, n=1,…,N by counting failures. However, we cannot obtain exactly Xn due to noises. In this case, we must establish new estimators for a and b only based on observations of M1,…,MN. In a simple heuristic, we may replace Xi by Mi in (8) to (9) and obtain estimators a* and b* for a and b, respectively. Based on simulation results in Section 5, we will conclude that a* and b* are bad estimators. In the following, we construct a good approximation of Xn based on M1,…,MN and obtain better estimators for a and b.Theorem 1.
M
=
{
M
n
:
n
=
1,2
,
…
} is a Poisson-faults model with an intensity X={Xn:n=1,2,…}. As N tends to ∞, with a probability 1, one has that
(10)1N-1∑n=1N-1‍(Mn-Xn)⟶0,1N-1∑n=1N-1‍(Mn+1Mn-Xn+1Xn)⟶0,1N-1∑n=1N-1‍(Mn2-Mn-Xn2)⟶0.See AppendixA for the proof of Theorem 1.Then, it follows from Theorem1 and (8) that we can solve linear equations
(11)A~(a~b~)=β~,A~=(∑n=1N-1‍Mn2-Mn,∑n=1N-1‍Mn∑n=1N-1‍Mn,N-1),β~=(∑n=1N-1‍Mn+1Mn∑n=1N-1‍Mn+1)
for estimators a~ and b~. Moreover, it follows from Theorem 1 and the convergence of a^ and b^ that a~ and b~ converge to a and b with probability 1, respectively. Therefore, as N tends to ∞, a~ and b~ converge to a and b with a probability 1, respectively. Similarly, from (9), we can construct an estimator
(12)σ~2=b~2+1N-1∑n=1N-1‍a~2(Mn2-Mn)+(Mn+12-Mn+1)-2a~MnMn+1-2b~Mn+1+2a~b~Mn
for the parameter σ2.


## 4. Bayes Estimate of the Instant Failure Rate

In this section, from observations ofM1,…,MN, we try to estimate the instant failure rate Xn. According to the general theory of Poisson compound decision problem [19], we establish a formula to the estimate. After investigating the time complexity of the formula, we turn to another quick algorithm.


### 4.1. Formula Involving Moments

Theorem 2.
M
=
{
M
1
,
…
,
M
N
} is a Poisson-fault model with the intensity X. (m1,…,mN) is a dataset of observations of M. Assume that the prior normal distribution of (X1,…,XN) has a mean μ~ and a covariance B. Let U=(U1,…,UN) be a normal random variable with mean μ=μ~-Be and covariance B, where e denotes N-dimensional vector whose every entry is 1. One has that, for n=1,…,N,
(13)E[Xn∣M1=m1,…,MN=mN]=E(Un∏i=1n‍Uimi)E∏i=1n‍Uimi.See AppendixB for the proof of Theorem 2.From Theorem2, we can obtain the Bayes estimate E[Xn∣M1,…,MN] by computing L+1th and Lth moments of a normal random variable U, where L=m1+⋯+mN is the total number of failures. The moment generating function ϕ(t) of U provides a tool to calculate the moments; see [20] for details of a moment generating function. For t=(t1,…,tn), the moment generating function ϕ(t) of U is defined by
(14)ϕ(t)=Eexp⁡(∑i=1N‍tiUi)=exp⁡(∑i=1N‍tiμi+12∑i=1N‍∑j=1N‍titjbi,j),
where (μ1,…,μN)=μ and (bi,j)=B. In general, nth derivative of ϕ(t) evaluated at t=0 equals a corresponding nth moment of U. Concretely, for nonnegative integers n1,…,nn,
(15)E(∏i=1N‍Uini)=∂(n1+⋯+nN)∂t1n1⋯∂tNnNϕ(t)|t=(0,…,0).On the other hand, differentiation ofϕ(t) with respect to tk yields
(16)∂∂tkϕ(t)=(μk+∑j=1N‍tjbk,j)ϕ(t).
Hence, for nonnegative integers s1,…,sn,
(17)∂(s1+⋯+sN)∂t1s1⋯∂tNsN(∂∂tkϕ(t))=(μk+∑j=1N‍tjbk,j)∂(s1+⋯+sN)∂t1s1⋯∂tNsNϕ(t)+∑j=1N‍sjbk,j∂(s1+⋯+sN-1)∂t1s1⋯∂tksk-1⋯∂tNsNϕ(t).
Moreover, it follows from (15) that
(18)E(Uk∏i=1N‍Uisi)=μkE(∏i=1N‍Uisi)+∑j=1N‍sjbk,jE(Ujsj-1∏i≠j‍Uisi).As we can obtain the1st and 2nd moments of U from its mean μ and covariance matrix B, the preceding (18) provides us with a formula to calculate any higher-order moment of U. However, the formula involves N order tensor. The time complexity of a corresponding algorithm for E(∏i=1N‍Uimi) is O(NL), where L=m1+⋯+mN.As an alternate way to calculateE(∏i=1N‍Uimi), we had studied a Monte Carlo method [15, 21]. In general, to achieve the specified accuracy of E(∏i=1N‍Uimi), the runtime of a Monte Carlo algorithm is proportional to the variance of the random variable ∏i=1N‍Uimi. Our simulation shows that the performance of the Monte Carlo algorithm is not good when L is slightly large.


### 4.2. Methodology of a Quick Algorithm

So we turn to another technique to calculate the Bayes estimate. The methodology is similar to [13, 14].Assume thatf(x1) and g(ϵn), n=2,3,…, are probability density functions (pdf) of random variables X1 and ϵn, n=2,3,…, respectively. As ϵn is independent with respect to ℱn-1, we have that the multivariate pdf of (X1,ϵ2,…,ϵN) is
(19)hN(x1,ϵ2,…,ϵN)=f(x1)∏n=2N‍g(ϵn).Suppose thatXn=xn, n=1,…,N, is a dataset of the intensity process. It follows from the definition of the Poisson-fault model that the conditional probability function of M1,…,MN is
(20)P(M1=m1,…,MN=mN∣X1=x1,…,XN=xN)=∏n=1N‍(xnmnmn!exp⁡(-xn)).
Therefore,
(21)P(M1=m1,…,MN=mN)=∫RN‍∏n=1N‍(xnmnmn!exp⁡(-xn))hNdτN.
Here hN=hN(x1,ϵ2,…,ϵN), dτN=dx1dϵ2⋯dϵN, and xn=axn-1+b+ϵn for n=2,…,N. From the above discussion, we can derive the following expectation of Xj conditional on a dataset of M1,…,MN(22)E[Xj∣M1=m1,…,MN=mN]=∫RN‍xj∏n=1N‍(xnmnmn!exp⁡(-xn))hNdτN.To compute integrals in (21) and (22), we study the following auxiliary function In(k1,…,kn) for n=1,…,N:
(23)In=∫Rn‍∏i=1n‍xikiexp⁡(-∑i=1n-1‍xi-1-aN-n+11-axn)hndτn.
As xn=axn-1+b+ϵn, writing cn=(1-bN-n+1)/(1-b), we have that
(24)In=exp⁡(-cnb)×∫Rn-1‍exp⁡(-cnaxn-1)hhhh×∏i=1n-1‍xikiexp⁡(-∑i=1n-1‍xi)hn-1hhhh×(∫R‍g(ϵn)(axn-1+b+ϵn)knexp⁡(-cnϵn)dϵn)hhhh×dτn-1.
Hence, an evolutionary equation of In follows from the binomial formula
(25)In(k1,…,kn)=e-cnb∑l=0kn‍CknlE[(b+ϵn)kn-le-cnϵn]hhhhhh×akn-lIn-1(k1,…,kn-1+l).That is, based on expectations of several simple functions ofx1 or ϵn, n=2,3,…, we can obtain iteratively In(k1,…,kn) from (25). Moreover, it follows from (21) and (22) that
(26)E[Xj∣M1=m1,…,MN=mN]=IN(m1,…,mj+1,…,mN)IN(m1,…,mN).


### 4.3. Algorithm

Comparing (26) with (18), we note that there is still a tensor IN(m1,…,mN) in (26). We try to simplify IN(m1,…,mN).We know thatϵn has a normal distribution N(0,σ2) and we assume that X1 has a prior normal distribution N(x1,σ12). Then we can derive a straightforward procedure to compute E[Xn∣M1=m1,…,MN=mN].Assume that a random variableY has a normal distribution N(μy,σy2). Following the method of Section 4.1, for any nonnegative integer k, we have
(27)E[(b+Y)ke-cY]=exp⁡(c2σy22-c)∑i=0k‍Cki(b+μy-cσy2)k-iJi.
Here J2i+1=0 and J2i=∏j=1i‍(2j-1) for i=0,1,….In the following, the right side of (27) is denoted by E(k,c,b,μy,σy2). Based on the method of Section 4.2, we can derive the following algorithm to compute I[k1,…,kN].To calculateIN(m1,…,mN), it is obvious that the time complexity of the preceding algorithm is O(NL2), where L=m1+⋯+mN. Moreover, we can obtain E[Xj∣M1=m1,…,MN=mN] from (26).


## 4.1. Formula Involving Moments

Theorem 2.
M
=
{
M
1
,
…
,
M
N
} is a Poisson-fault model with the intensity X. (m1,…,mN) is a dataset of observations of M. Assume that the prior normal distribution of (X1,…,XN) has a mean μ~ and a covariance B. Let U=(U1,…,UN) be a normal random variable with mean μ=μ~-Be and covariance B, where e denotes N-dimensional vector whose every entry is 1. One has that, for n=1,…,N,
(13)E[Xn∣M1=m1,…,MN=mN]=E(Un∏i=1n‍Uimi)E∏i=1n‍Uimi.See AppendixB for the proof of Theorem 2.From Theorem2, we can obtain the Bayes estimate E[Xn∣M1,…,MN] by computing L+1th and Lth moments of a normal random variable U, where L=m1+⋯+mN is the total number of failures. The moment generating function ϕ(t) of U provides a tool to calculate the moments; see [20] for details of a moment generating function. For t=(t1,…,tn), the moment generating function ϕ(t) of U is defined by
(14)ϕ(t)=Eexp⁡(∑i=1N‍tiUi)=exp⁡(∑i=1N‍tiμi+12∑i=1N‍∑j=1N‍titjbi,j),
where (μ1,…,μN)=μ and (bi,j)=B. In general, nth derivative of ϕ(t) evaluated at t=0 equals a corresponding nth moment of U. Concretely, for nonnegative integers n1,…,nn,
(15)E(∏i=1N‍Uini)=∂(n1+⋯+nN)∂t1n1⋯∂tNnNϕ(t)|t=(0,…,0).On the other hand, differentiation ofϕ(t) with respect to tk yields
(16)∂∂tkϕ(t)=(μk+∑j=1N‍tjbk,j)ϕ(t).
Hence, for nonnegative integers s1,…,sn,
(17)∂(s1+⋯+sN)∂t1s1⋯∂tNsN(∂∂tkϕ(t))=(μk+∑j=1N‍tjbk,j)∂(s1+⋯+sN)∂t1s1⋯∂tNsNϕ(t)+∑j=1N‍sjbk,j∂(s1+⋯+sN-1)∂t1s1⋯∂tksk-1⋯∂tNsNϕ(t).
Moreover, it follows from (15) that
(18)E(Uk∏i=1N‍Uisi)=μkE(∏i=1N‍Uisi)+∑j=1N‍sjbk,jE(Ujsj-1∏i≠j‍Uisi).As we can obtain the1st and 2nd moments of U from its mean μ and covariance matrix B, the preceding (18) provides us with a formula to calculate any higher-order moment of U. However, the formula involves N order tensor. The time complexity of a corresponding algorithm for E(∏i=1N‍Uimi) is O(NL), where L=m1+⋯+mN.As an alternate way to calculateE(∏i=1N‍Uimi), we had studied a Monte Carlo method [15, 21]. In general, to achieve the specified accuracy of E(∏i=1N‍Uimi), the runtime of a Monte Carlo algorithm is proportional to the variance of the random variable ∏i=1N‍Uimi. Our simulation shows that the performance of the Monte Carlo algorithm is not good when L is slightly large.


## 4.2. Methodology of a Quick Algorithm

So we turn to another technique to calculate the Bayes estimate. The methodology is similar to [13, 14].Assume thatf(x1) and g(ϵn), n=2,3,…, are probability density functions (pdf) of random variables X1 and ϵn, n=2,3,…, respectively. As ϵn is independent with respect to ℱn-1, we have that the multivariate pdf of (X1,ϵ2,…,ϵN) is
(19)hN(x1,ϵ2,…,ϵN)=f(x1)∏n=2N‍g(ϵn).Suppose thatXn=xn, n=1,…,N, is a dataset of the intensity process. It follows from the definition of the Poisson-fault model that the conditional probability function of M1,…,MN is
(20)P(M1=m1,…,MN=mN∣X1=x1,…,XN=xN)=∏n=1N‍(xnmnmn!exp⁡(-xn)).
Therefore,
(21)P(M1=m1,…,MN=mN)=∫RN‍∏n=1N‍(xnmnmn!exp⁡(-xn))hNdτN.
Here hN=hN(x1,ϵ2,…,ϵN), dτN=dx1dϵ2⋯dϵN, and xn=axn-1+b+ϵn for n=2,…,N. From the above discussion, we can derive the following expectation of Xj conditional on a dataset of M1,…,MN(22)E[Xj∣M1=m1,…,MN=mN]=∫RN‍xj∏n=1N‍(xnmnmn!exp⁡(-xn))hNdτN.To compute integrals in (21) and (22), we study the following auxiliary function In(k1,…,kn) for n=1,…,N:
(23)In=∫Rn‍∏i=1n‍xikiexp⁡(-∑i=1n-1‍xi-1-aN-n+11-axn)hndτn.
As xn=axn-1+b+ϵn, writing cn=(1-bN-n+1)/(1-b), we have that
(24)In=exp⁡(-cnb)×∫Rn-1‍exp⁡(-cnaxn-1)hhhh×∏i=1n-1‍xikiexp⁡(-∑i=1n-1‍xi)hn-1hhhh×(∫R‍g(ϵn)(axn-1+b+ϵn)knexp⁡(-cnϵn)dϵn)hhhh×dτn-1.
Hence, an evolutionary equation of In follows from the binomial formula
(25)In(k1,…,kn)=e-cnb∑l=0kn‍CknlE[(b+ϵn)kn-le-cnϵn]hhhhhh×akn-lIn-1(k1,…,kn-1+l).That is, based on expectations of several simple functions ofx1 or ϵn, n=2,3,…, we can obtain iteratively In(k1,…,kn) from (25). Moreover, it follows from (21) and (22) that
(26)E[Xj∣M1=m1,…,MN=mN]=IN(m1,…,mj+1,…,mN)IN(m1,…,mN).


## 4.3. Algorithm

Comparing (26) with (18), we note that there is still a tensor IN(m1,…,mN) in (26). We try to simplify IN(m1,…,mN).We know thatϵn has a normal distribution N(0,σ2) and we assume that X1 has a prior normal distribution N(x1,σ12). Then we can derive a straightforward procedure to compute E[Xn∣M1=m1,…,MN=mN].Assume that a random variableY has a normal distribution N(μy,σy2). Following the method of Section 4.1, for any nonnegative integer k, we have
(27)E[(b+Y)ke-cY]=exp⁡(c2σy22-c)∑i=0k‍Cki(b+μy-cσy2)k-iJi.
Here J2i+1=0 and J2i=∏j=1i‍(2j-1) for i=0,1,….In the following, the right side of (27) is denoted by E(k,c,b,μy,σy2). Based on the method of Section 4.2, we can derive the following algorithm to compute I[k1,…,kN].To calculateIN(m1,…,mN), it is obvious that the time complexity of the preceding algorithm is O(NL2), where L=m1+⋯+mN. Moreover, we can obtain E[Xj∣M1=m1,…,MN=mN] from (26).


## 5. Simulation Results

We will study the performance of a power transformer from the viewpoint of a Poisson-fault model. It is reckless to apply our method to an empirical test: it is difficult to explore all features of a power transformer based on empirical datasets due to the loss of the intensityXn, and hence we cannot evaluate the performance of the method. Through simulations, we may obtain observations of both Xn and Mn. In the section, we focus on simulation experiments about failures of a power transformer.


### 5.1. Dataset

The dataset should be as similar as possible to real world dataset. In the following, we start from a model used in [12] for transformers to deduce reasonable a, b, ϵn, and X1 for (3). In the model used in [12], the relative average lifetime
(28)t=A22-(T-98)/6=82371.3A2exp⁡(-0.1155T).
Here T is the hottest spot temperature of a power transformer in degrees Celsius. A2 is the average lifetime of a transformer at T=98. A typical A2 is two decades. Comparing (28) to (6), we have that k=0.1155. It follows from (7) that
(29)a=e-kΔ.
Here the benchmark interval Δ is a change in temperature between two successive working conditions. In the following simulations, we choose Δ=3 and a=0.707. Besides, the level n is from 1 to 9, corresponding to the temperature T=134,…,110, respectively. That is, the total number of levels is N=9.Now we turn to the choices ofb and ϵi. b represents the bedrock errors which does not appear in the Arrhenius model. As the Arrhenius model is a widely accepted model, b is a small number. In the following simulations, we use different settings 0.01, 0.02, 0.05, and 0.1 for b. ϵn represents a difference of the behavior between a real transformer and the Arrhenius model. It is also a small random variable. In following simulations, we assume that ϵn has a normal distribution N(0,0.01).X
1 is the average of the number of failures at level 1, which corresponds to the most severe working condition. We assume that T=134 at level 1. Then the average lifetime t=0.0156A2 from (28). For a typical A2=20, it follows from (5) that
(30)X1=dp0.312.dp may have a wide range for choosing. In the following simulations, we set dp=8. Besides, the model allows a small noise for X1. Hence, we assume that X1 has a normal distribution N(25.64,0.01).The simulation ofMn based on Xn is according to the classical method to simulation of a Poisson process; see [13] for example.


### 5.2. Simulation Results

Based on datasets provided by the method introduced in Section5.1, we perform three experiments to test performances of the model and the method introduced in the above two sections.In the first experiment, we compare estimatorsa~ and b~ with a* and b*. For different value of b, we repeat simulations 400 times to obtain 400 datasets. We present the means and standard errors of a~ and b~ of these datasets in Table 1. We also present the means and standard errors of a* and b* in Table 2. Noting that the true value of a is 0.707 from (29), Tables 1 and 2 show that the standard errors of both a~ and a* are small, but the mean of a~ is closer to the true value than that of a*. Moreover, from these tables, we may conclude that the performances of both the mean and the standard errors of b~ are better than that of b*.Table 1
Mean and standard errors ofa~ and b~.
True value ofb
Mean ofa~
Standard errors ofa~
Mean ofb~
Standard errors ofb~
0.01
0.7039
0.2025
0.0308
0.0711
0.02
0.7202
0.1881
0.0197
0.0806
0.05
0.7137
0.1996
0.0465
0.0852
0.1
0.7024
0.1971
0.1075
0.0738Table 2
Mean and standard errors ofa* and b*.
True value ofb
Mean ofa*
Standard errors ofa*
Mean ofb*
Standard errors ofb*
0.01
0.6067
0.1384
1.0252
1.1758
0.02
0.5941
0.1374
1.1161
1.1550
0.05
0.6062
0.1371
1.0964
1.1519
0.1
0.5969
0.1446
1.2584
1.2881In the second experiment, we setb=0.01 and run 400 simulations. For every dataset, based on the estimators a~ and b~ for the dataset, we deduce by our algorithm the expectations E[Xn∣M1,…,M9] for n=1,…,9. In Figure 1, we compare these expectations with Xn, n=1,…,9. As there are 400 datasets of Xn, the benchmark curve is the mean of Xn. We may conclude from Figure 1 that the mean of expectations E[Xn∣M1,…,M9] is close to the true value and the standard errors are small.Figure 1
The mean and the standard errors of expectations. The solid curve is the benchmark curve. The dashed curves stand for the mean ofE[Xn∣M1,…,M9] with standard errors. And the stars are the mean of E[Xn∣M1,…,M9].The third experiment concerns the time complexity of the method. We randomly selected12 groups of nonnegative m1,…,m9 such that L=m1+…+m9 is 3,4,… or 14. We use three methods to calculate E[X5∣M1=m1,…,M9=m9]: the algorithm in Section 4.3, formula (18), and a Monte Carlo simulation suggested in Section 4.1. The Monte Carlo simulation stops when runtime exceeds half an hour or the standard error is less than 0.1. The outputs presented in Figures 2 and 3 show the corresponding runtime.Figure 2
The Bayes estimateE[X5∣M1=m1,…,M9=m9]. The solid curve is from (18). The dashed curve with circles stands for the algorithm in Section 4.3. The dashed curve with crosses is the output of the Monte Carlo simulation.Figure 3
Runtime for three methods. The solid curve with dots stands for (18). The dashed curve with circles stands for the algorithm in Section 4.3. The dashed curve with crosses stands for the Monte Carlo simulation.For smallL, the result of Bayes estimate are similar for three methods. The performance of the Monte Carlo simulation is not good for L≥11 because its runtime exceeds its limit. On the other hand, the runtime of formula (18) is unacceptable for L=9. In fact, we do not test (18) for L≥10 due to its long runtime at L=9.


## 5.1. Dataset

The dataset should be as similar as possible to real world dataset. In the following, we start from a model used in [12] for transformers to deduce reasonable a, b, ϵn, and X1 for (3). In the model used in [12], the relative average lifetime
(28)t=A22-(T-98)/6=82371.3A2exp⁡(-0.1155T).
Here T is the hottest spot temperature of a power transformer in degrees Celsius. A2 is the average lifetime of a transformer at T=98. A typical A2 is two decades. Comparing (28) to (6), we have that k=0.1155. It follows from (7) that
(29)a=e-kΔ.
Here the benchmark interval Δ is a change in temperature between two successive working conditions. In the following simulations, we choose Δ=3 and a=0.707. Besides, the level n is from 1 to 9, corresponding to the temperature T=134,…,110, respectively. That is, the total number of levels is N=9.Now we turn to the choices ofb and ϵi. b represents the bedrock errors which does not appear in the Arrhenius model. As the Arrhenius model is a widely accepted model, b is a small number. In the following simulations, we use different settings 0.01, 0.02, 0.05, and 0.1 for b. ϵn represents a difference of the behavior between a real transformer and the Arrhenius model. It is also a small random variable. In following simulations, we assume that ϵn has a normal distribution N(0,0.01).X
1 is the average of the number of failures at level 1, which corresponds to the most severe working condition. We assume that T=134 at level 1. Then the average lifetime t=0.0156A2 from (28). For a typical A2=20, it follows from (5) that
(30)X1=dp0.312.dp may have a wide range for choosing. In the following simulations, we set dp=8. Besides, the model allows a small noise for X1. Hence, we assume that X1 has a normal distribution N(25.64,0.01).The simulation ofMn based on Xn is according to the classical method to simulation of a Poisson process; see [13] for example.


## 5.2. Simulation Results

Based on datasets provided by the method introduced in Section5.1, we perform three experiments to test performances of the model and the method introduced in the above two sections.In the first experiment, we compare estimatorsa~ and b~ with a* and b*. For different value of b, we repeat simulations 400 times to obtain 400 datasets. We present the means and standard errors of a~ and b~ of these datasets in Table 1. We also present the means and standard errors of a* and b* in Table 2. Noting that the true value of a is 0.707 from (29), Tables 1 and 2 show that the standard errors of both a~ and a* are small, but the mean of a~ is closer to the true value than that of a*. Moreover, from these tables, we may conclude that the performances of both the mean and the standard errors of b~ are better than that of b*.Table 1
Mean and standard errors ofa~ and b~.
True value ofb
Mean ofa~
Standard errors ofa~
Mean ofb~
Standard errors ofb~
0.01
0.7039
0.2025
0.0308
0.0711
0.02
0.7202
0.1881
0.0197
0.0806
0.05
0.7137
0.1996
0.0465
0.0852
0.1
0.7024
0.1971
0.1075
0.0738Table 2
Mean and standard errors ofa* and b*.
True value ofb
Mean ofa*
Standard errors ofa*
Mean ofb*
Standard errors ofb*
0.01
0.6067
0.1384
1.0252
1.1758
0.02
0.5941
0.1374
1.1161
1.1550
0.05
0.6062
0.1371
1.0964
1.1519
0.1
0.5969
0.1446
1.2584
1.2881In the second experiment, we setb=0.01 and run 400 simulations. For every dataset, based on the estimators a~ and b~ for the dataset, we deduce by our algorithm the expectations E[Xn∣M1,…,M9] for n=1,…,9. In Figure 1, we compare these expectations with Xn, n=1,…,9. As there are 400 datasets of Xn, the benchmark curve is the mean of Xn. We may conclude from Figure 1 that the mean of expectations E[Xn∣M1,…,M9] is close to the true value and the standard errors are small.Figure 1
The mean and the standard errors of expectations. The solid curve is the benchmark curve. The dashed curves stand for the mean ofE[Xn∣M1,…,M9] with standard errors. And the stars are the mean of E[Xn∣M1,…,M9].The third experiment concerns the time complexity of the method. We randomly selected12 groups of nonnegative m1,…,m9 such that L=m1+…+m9 is 3,4,… or 14. We use three methods to calculate E[X5∣M1=m1,…,M9=m9]: the algorithm in Section 4.3, formula (18), and a Monte Carlo simulation suggested in Section 4.1. The Monte Carlo simulation stops when runtime exceeds half an hour or the standard error is less than 0.1. The outputs presented in Figures 2 and 3 show the corresponding runtime.Figure 2
The Bayes estimateE[X5∣M1=m1,…,M9=m9]. The solid curve is from (18). The dashed curve with circles stands for the algorithm in Section 4.3. The dashed curve with crosses is the output of the Monte Carlo simulation.Figure 3
Runtime for three methods. The solid curve with dots stands for (18). The dashed curve with circles stands for the algorithm in Section 4.3. The dashed curve with crosses stands for the Monte Carlo simulation.For smallL, the result of Bayes estimate are similar for three methods. The performance of the Monte Carlo simulation is not good for L≥11 because its runtime exceeds its limit. On the other hand, the runtime of formula (18) is unacceptable for L=9. In fact, we do not test (18) for L≥10 due to its long runtime at L=9.


## 6. Conclusions

In this paper, we proposed a Poisson-fault model to describe the failure of a power transformer at different working conditions. The model is an extension to the Arrhenius model. We modified the traditional time series techniques and proposed a framework to calculate parameters of the model. After exploring three methods with different time complexity, we established a quick algorithm to estimate the failure rate of a power transformer. The algorithm is tested on simulation datasets for a power transformer. By simulation studies, we have shown that the estimators we proposed have better performance than the estimators based on the classical method of time series. Moreover, the runtime of our method is the best among three methods. The estimator of the intensity is near the true value and has a small standard error.From the viewpoint of an engineer, we use a slightly ideal assumption that the change of two successive working conditions is constant. If the change is not a constant, we do not know whether the runtime of the proposed method is still acceptable to applications. On the other hand, from the viewpoint of statistician, it is an unnatural assumption that a parameter of the model has a prior normal distribution, because the posterior distribution is not a normal distribution. To extend the proposed method to a case with general change and to establish a natural prior distribution for a parameter are future research directions. Ongoing work is investigating.



---

*Source: 945258-2014-05-14.xml*