Linear Arrhenius-Weibull Model for Power Transformer Thermal Stress Assessment

Arrhenius equations with Weibull distribution have been broadly deployed to quantify the loss of life and probability of failure of power transformers. This model is highly nonlinear, and this non-linearity makes it challenging to use this model for grid management and optimization. In this study, equations are linearized using Taylor series expansion to provide a linear formulation for transformer loss of life and probability of failure as a function of ambient temperature and transformer loading. The proposed formulation allows transformer constraints to be incorporated in grid management applications within conventional optimization approaches, such as linear programming, and decreases the calculation burden caused by nonlinearity.


I. INTRODUCTION
Electric power distribution transformers are significant elements of power systems, and their operation plays an important role in the reliability and resilience of the system. Hence, estimating the stress caused by their operating conditions is important for mitigating their loss of life and the probability of failure.
Life of the power transformers mainly depends on the life of their cellulose paper insulation which deteriorates with thermal stress, moisture, and oxygen. Oxygen and moisture can be managed to remain minimal; however, temperature is not as manageable [1]. The life of the transformer is governed by its thermal behavior [2]. Because the temperature in a transformer is distributed non-linearly, the insulation exposed to the hottest spot temperature experiences the fastest degradation and would be the most likely to fail [3]. Hence, using the hottest spot temperature in thermal models provides the best estimation for transformer life.

A. MOTIVATION
Thermal stress is one of the main reasons for transformers' loss of life and failure [4]. The temperature of transformer windings and insulation systems, oil for example, is a major The associate editor coordinating the review of this manuscript and approving it for publication was Arpan Kumar Pradhan.
factor effecting transformer loading capability [5]. Quantifying transformers' loss of life and probability of failure as a function of loading may be then used for grid planning and operation studies.
For planning purposes, the calculation burden may not be as important, but for operational applications, it is important to reduce the calculation burden as much as possible due to the limited time for assessing the situation and making an appropriate decision.
For both planning and operational purposes, transformer thermal stress assessment may be utilized in solving an optimization problem. There are simple and effective algorithms for solving large scale linear optimization problems using a low-cost processor. However, if the formulation is non-linear, finding an optimal solution may become a challenging and computationally expensive task. Providing a linear formulation for transformer loss of life and probability of failure enables the solution of planning and operational optimization problems using standard approaches and processors.

B. PRIOR RESEARCH
Quantifying the loss of life and the probability of failure for transformers has been a subject of extensive research. In one of the first publications in this area, historical failure data were used to calculate the probability distribution function VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ using statistical techniques [6]. Because sufficient failure data for a transformer may not be available and the lifetime of the transformer is relatively long, studies have attempted to provide a more accurate estimation of the parameters of the distribution functions [7], [8]. Several distribution functions have been studied for the purpose of transformer life evaluation. Some have used normal and log-normal distributions [7], [9], while another used Perk's hazard function [10]. Yet another paper used generalized exponential distribution [11]. Weibull distribution has been discussed extensively for this purpose [7], [9], [12]- [16]. All distributions used in these methods are highly nonlinear. Arrhenius model is used for life assessment of various power system components such as transformers [17], surge arresters [18], [19], electrical insulating materials [20], and cables [21].
Temperature related chemical changes in the electrical insulation were discussed in [22] for the first time. The life of an insulator can be impacted by the surrounding insulation medium as well [23]. The surrounding medium for the oil-filled transformers is oil and thermal stress may cause chemical reactions that lead to release of several gases in the transformer oil [24]. These gases can be analyzed by chromatography and used to detect thermal stress [25]. Various ratio-based methods are introduced in several standards including IEEE C57.104-2008 [26], IEC 60599 [27], and ASTM D3612-90 [28] as well as graphical methods [29], [30]. In addition to the conventional methods, there have been recent efforts for dissolved gas analysis as described in [31] and [32].

C. CONTRIBUTION OF THIS PAPER
We adopted the formulation of the Arrhenius model for life assessment with a Weibull distribution. The linear form of this model was analytically obtained using Taylor series; the proposed linear form is used for life evaluation in a Monte Carlo simulation of 250,000 cases with different loading scenarios. We show that the proposed formulation has acceptable performance compared to the original nonlinear form.
The remainder of this paper is organized as follows. The formula for calculating the hottest spot temperature of the transformers is explained in Section II. In Section III, the Arrhenius model is discussed. Linearization of the loss of life and probability of failure equations are described in Sections IV and V, respectively. The validation and conclusions are outlined in Sections VI and VII, respectively, followed by references.

II. HOTTEST SPOT TEMPERATURE
The hottest spot temperature is the main factor in assessing the loss of life and the probability of failure of a transformer. The deterioration of transformer insulation is a function of time and temperature [5], [33]. There are several IEEE standards for calculating the hottest spot temperature [34]. The hottest spot temperature, θ H , consists of three parts: the ambient temperature (θ a ) in Kelvin, the top oil temperature rise over ambient temperature ( θ TO ), and the conductor temperature rise over the top oil temperature ( θ W ). It is calculated as follows: Top oil temperature rise over ambient temperature can be calculated using (2): where θ TO,rated is the transformer's top oil temperature rise over the ambient temperature under nominal loading, x is the ratio of the apparent loading of the transformer to its nameplate rating, R is the ratio of loss at rated load to no load loss, and n is the exponential power of loss versus the top oil temperature rise. We rewrite (2) as where The conductor temperature rise over top oil temperature can be calculated using (6): where θ W ,rated is the conductor temperature rise over the top oil temperature at the rated load, and m is the exponential power of the winding loss versus winding gradient. The equations introduced so far are used for steady-state calculations. To calculate the transient temperature, the transient values of the top oil temperature rise over the ambient temperature and the conductor temperature rise over the top oil temperature can be calculated using (7) and (8), respectively: where ini and ult refer to the steady-state values at the beginning and end of each time step, respectively. Terms τ W and τ TO are the conductor and oil thermal time constants, respectively, and t is the time passed from the beginning of the time step. In our study, it is assumed that the time steps are long enough for the temperatures to reach steady state values at the end of each time step. Thus, the initial values of each time step are the ultimate values of the previous time step. A unit of time can be defined for the time interval of each temperature evaluation, based on the loading of the transformer. Whatever interval is assumed as the unit of time, the oil and conductor thermal time constants should be adjusted accordingly. To simplify the problem and to avoid time as a variable in the hottest spot temperature, the average temperature rise over a unit time interval is calculated using (9) and (10): where T is the time at the start of the time step. Solving the integrals in (9) and (10), the temperature rises can be obtained using (11) and (12): Substituting the initial and ultimate temperatures using (3) and (6) results in what is shown in (13) and (14): where Using (1), (13), and (14), the hottest spot temperature is a function of the loading in the present time step and the previous time step, as shown in (17).

III. ARRHENIUS MODEL
The Arrhenius model is a well-known life-stress model when the acceleration variable is thermal. The Arrhenius reaction rate equation is where R is the reaction rate; that is, the speed of reaction, C is a nonthermal constant, E a is the activation energy (eV), k B is Boltzmann's constant (eV/K), and θ is the absolute temperature (K). Life is proportional to the inverse action rate; thus, λ can be defined as a quantifiable life measure using (18), as shown in (19): where A and B are empirical parameters. For the electric transformer age assessment, the hottest spot temperature is considered as the temperature in the Arrhenius model [34]. The Arrhenius equation for the transformer life-stress model is shown in (20):

IV. TRANSFORMER LOSS OF LIFE
The model in Equation (20) can be used along with its value for normal operating conditions to determine the per-unit life that the transformer loses in the time step. This is called the aging acceleration factor (F AA ), as shown in (21): where θ H ,ref is the reference temperature at which the transformer is supposed to operate normally. To simplify (21), the normal aging part is replaced with a constant D in (22): Now, we attempt to find a linear approximation for the acceleration factor. The chain rule in Leibniz's notation is given in (23): The two-variable Taylor series expansion for f ( (24): By removing the higher-order terms (h.o.t.) in equation (24), the function can be estimated as shown in (25): where x * 1 and x * 2 are the initial and ultimate loading conditions for the transformer, respectively. Using (22), the derivative of the acceleration factor with respect to x 1 VOLUME 10, 2022 and x 2 can be calculated as shown in (26) and (27), respectively: Using (22), the derivative of the acceleration factor with respect to the hottest spot temperature can be calculated as shown in (28): Using (13), the derivatives of the top oil temperature rise over the ambient temperature with respect to x 1 and x 2 can be calculated as shown in (29) and (30): Using (14), the derivatives of the conductor temperature rise over the top oil temperature with respect to x 1 and x 2 can be calculated as shown in (31) and (32): Using (17), (29), and (31), the derivative of the hottest spot temperature with respect to x 1 is calculated as shown in (33): Similarly, using (17), (30), and (32), the derivative of the hottest spot temperature with respect to x 2 is calculated as shown in (34): Using (26), (28), and (33), the derivative of the acceleration factor with respect to x 1 is calculated as shown in (35): Similarly, using (27), (28), and (33), the derivative of the acceleration factor with respect to x 2 is calculated as shown in (36): Using (25), (35), and (36), and simplifying the notation θ H (x * 1 , x * 2 ) to θ H , the acceleration factor can be estimated in linear form for x 1 and x 2 as shown in (37): The remaining question to be addressed is the values of x * 1 and x * 2 . Equation (37) is a general equation and based on the study and different operating conditions, it may be more effective to use the predicted loading values for x * 1 and x * 2 . In addition, for situations in which the loading may vary significantly, several values of x * 1 and x * 2 can be used in order to construct a linear approximation. Hence, simulations described in Section IV explore multiple choices of values of x * 1 and x * 2 in (37).

V. TRANSFORMER PROBABILITY OF FAILURE
It is suggested in [5] that the Weibull distribution is suitable for use in the transformer's probability of failure in its life span. In addition, as discussed in Section III, the Arrhenius model is used to model the thermal degradation of the transformer's insulation. Thus, the probability of a thermal-related failure for a transformer can be quantified using Arrhenius-Weibull model. The probability density function (PDF) for two-parameter Weibull distribution is given by (38): where δ and λ are the shape and characteristic life parameters, respectively. As explained in [5], [14], and [16], the characteristic life parameter is the same as the Arrhenius relationship given by (20). The cumulative distribution function (CDF) of a continuous random variable X can be calculated from its PDF f X (t) as shown in (39): Using (38) and (39), the CDF of the Arrhenius-Weibull model is as shown in (40): (40) By substituting λ from (20) in (40), the CDF can be written in the form shown in (41): (41) Reference [35] proposes a formulation for the probability of failure, as shown in (42): where t is the length of the time interval, and f (t, x 1 , x 2 ) is the failure probability density function. Using (42), the failure probability can be written as shown in (43).
Using (38) and (40), equation (43) can be written as shown in (44): (44) To linearize (44), the method explained in Section IV is deployed. Using (23), the derivative of the transformer probability of failure with respect to x 1 and x 2 can be calculated as shown in (45) and (46), respectively: Using (44), the derivative of the transformer probability of failure with respect to λ can be calculated as shown in (47): Using (20), the derivative of λ with respect to the hottest spot temperature can be calculated as shown in (48): The derivative of the hottest spot temperature with respect to x 1 and x 2 are calculated in Section IV and shown in (33) and (34). Using (33), (47), and (48), the derivative of the transformer probability of failure with respect to x 1 is shown in (49): (49) Using (34), (47), and (48), the derivative of the transformer probability of failure with respect to x 1 is shown in (50): (50) Using (25), (49), and (50), the transformer probability of failure can be estimated in linear form for x 1 and x 2 as shown in (51): The values of x * 1 and x * 2 can be assigned similar to the explanation in this regard in the last paragraph of Section IV. Considering that failure occurs only once, the probability of failure is conditional on the fact that the transformer has survived until the time of analysis [16]. This type of probability depends on the behavior of the equipment in the past; therefore, the Markov model is inapplicable. This formulation is based on a fixed operating temperature. Hence, to deploy it to evaluate the probability of failure for each moment, it should be calculated for similar aging conditions, that is, the equivalent operating time if the transformer had operated in the current operating condition from the beginning. For example, if the operating condition is normal, the transformer's effective aging in one year is equal to one year. However, if the transformer operates under the overloading condition, its effective aging may reach one year in ten days. Hence, this formulation should use the operating time for the operating condition that leads to the same aging condition that the transformer is in rather than the real operating time. The effective age of the transformer can be calculated using (20), as shown in (52): In (52), T is the equivalent age of the transformer if it was operating in the operating condition of the current time-step from when it was installed. T can be calculated as shown in (53): . (53) Using (20), equation (53) can be written as shown in (54): where T eq is defined as shown below:  It should be noted that this time conversion should only be applied when the stress level is greater than the stress at the normal level and the F AA is greater than one. Under this condition, equation (44) can be written as shown in (56): (56)  Using (56), the derivative of the transformer probability of failure with respect to λ is modified as shown in (57): Hence, the linearized probability of failure when F AA >1 is shown in (58): As can be seen in Tables 3 and 4, the relative error of deploying the proposed linear approximation is acceptable.

VI. VALIDATION
In this section, we describe a Monte Carlo simulation to evaluate the precision of the developed linear equations in comparison to the non-linear equations for transformer loss of life and probability of failure. For this purpose, five levels of loading for the transformer as shown in Table 1 were considered.
The values x * shown in Table 1 are the linearization points that will be used in (42) and (56) as x * 1 and x * 2 , which can be employed to obtain a linear approximation. Because five levels of loading are considered, based on the initial and ultimate loads, there are 25 linear pieces. Although in this paper all scenarios are considered, in the real world, it is not likely that the loading of a transformer jumps from ''Low'' to ''Extreme'', for example, so many of the 25 defined linear approximations may not be used. However, they are considered for the purpose of generality. The characteristics of the transformer used in this case study are adopted from [14] and are shown in Table 2.
In the Monte Carlo simulation, 10,000 scenarios of loading are generated for each of the 25 linear pieces. The sample  Table 1 are shown in blue.
points are generated randomly in the specified ranges for loading. In this simulation, the maximum limit for x 1 and x 2 values is 5/3 and the minimum limit is zero. The average relative error of linearizing in comparison to the non-linear Arrhenius-Weibull model for transformer acceleration factor and probability of failure for all 10,000 scenarios are calculated and shown in Table 3 and Table 4, respectively. The average relative errors of the linear acceleration factor and VOLUME 10, 2022 the probability of failure are calculated using (59) and (60): To visualize the accuracy of the proposed formulations for transformer loss of life and probability of failure, a new simulation is performed and x 2 is drawn with respect to x 1 for fixed values of F AA and P f . The results are shown in Figures 1 and 2. In these figures, the red lines represent the original non-linear models. To draw the red lines, the values of F AA and P f are assumed to be fixed and 160 uniformly distributed points from 0.01 to 1.6 are assigned to x 1 . For each sample points of x 1 and fixed value of F AA or P f , the value of x 2 is calculated using the non-linear Arrhenius-Weibull model. The 32 blue points in each figure are uniformly distributed for x 1 ∈ [0.05,1.6]. Similar to the red line, for each sample point of x 1 and fixed value of F AA or P f , the value of x 2 is calculated using the corresponding linear approximation proposed in this paper. The linear formulation is calculated based on (x * 1 ,x * 2 ) point representing the interval containing the corresponding red point, as shown in Table 1. As seen in the figures, the proposed linear model is close to the original nonlinear model.

VII. CONCLUSION
This paper proposes a linear form of the Arrhenius-Weibull model for transformer loss of life and the probability of failure. Linearization was performed using a Taylor series under different loading conditions. The proposed formulation is used in a case study that uses a Monte Carlo simulation of 250,000 cases. The average errors of the linear models compared to the non-linear model are calculated, and it is observed that the error is acceptable, which makes the proposed formulation accurate enough to be deployed. MLADEN KEZUNOVIC (Life Fellow, IEEE) is currently a Regents Professor at Texas A&M University. He is also a Principal of XpertPower Associates TM , a consulting firm specializing in power systems data analytics for over 30 years. His expertise is in protective relaying, computational intelligence for automated power system disturbance analysis, predictive data analytics for outage management, and smart grids. He authored over 650 articles and research reports, ten books and book chapters, gave over 150 keynotes and invited talks, presented over 50 seminars and short courses, and consulted for over 50 companies worldwide. He is a CIGRE Fellow, a Honorary, and a Distinguished Member. He is a Registered Professional Engineer in Texas.
SERGIY BUTENKO received the B.S. and M.S. degrees in mathematics from the Taras Shevchenko National University of Kyiv, Ukraine, in 1998 and 1999, respectively, and the M.S. and Ph.D. degrees in industrial engineering from the University of Florida, in 2001 and 2003, respectively. He is currently a Professor with the Wm. Michael Barnes '64 Department of Industrial and Systems Engineering, Texas A&M University. His research interests include fundamental and applied aspects of global and discrete optimization, as well as network analysis. Since 2013, he has been serving as the Editor-in-Chief for the Journal of Global Optimization. VOLUME 10, 2022