Estimating Wind Farm Transformers Rating through Lifetime Characterization Based on Stochastic Modeling of Wind Power

Abstract

This paper deals with the problem of the optimal rating of mineral-oil-immersed transformers in large wind farms. The optimal rating is derived based on the probabilistic analyses of wind power generation through the Ornstein–Uhlenbeck stochastic process and on thermal model of the transformer through the integration of stochastic differential equations. These analyses allow the stochastic characterization of lifetime reduction of the transformer and then its optimal rating through a simple closed form. The numerical application highlights the effectiveness and easy applicability of the proposed methodology. The proposed methodology allows deriving the rating of transformers which better fits the specific peculiarities of wind power generation. Compared to the conventional approaches, the proposed method can better adapt the transformer size to the intermittence and variability of the power generated by wind farms, thus overcoming the often-recognized reduced lifetime.

1. Introduction

Due to the large use of transformers, their accurate rating is crucial to avoid faults and guarantee secure and reliable operation of power systems. The task of transformer rating is typically performed through the estimation of the lifetime duration under a predicted loading condition. Lifetime duration strictly depends on the degradation of solid insulation, which is the primary reason why a transformer reaches the end of its life. Continuous degradation mainly occurs because of thermal stresses which affect transformers as a consequence of the heat generated during operation [1]. Thermal models are available in the literature and are used in consolidated international standards, such as the IEC Std 60076 [2]. Based on the loading current, they provide accurate estimations of the temperatures affecting the heat-sensitive parts of the transformers which are used in lifetime calculation.

Methods for transformer rating consist of choosing the size of the transformer which better fit the expected current flowing in the transformer with the desired lifetime duration, thereby avoiding accelerated degradation. The key aspect to optimize the determination of the transformer size relies in the determination of loading condition, which is a very complex task in modern power systems being subject to large uncertainties and random time disturbances, which affect all aspects of the operation of the system. Fluctuation in the power demand and intermittent nature of renewable generation are the most significant examples. It is also well known how challenging is to plan both the electric grid and the various components due to continuous power variations. In the case of transformers, aspects of uncertainties need to be carefully considered in order to accurately determine loading conditions [3].

This paper focuses on an accurate method for the optimal rating of oil-immersed power transformers connected to wind farms. The effect of fast variations of wind velocity implies that the transformer can be repeatedly overloaded without impacting on the loading current mean value but causing overheating, thus implying an erroneous sizing of the transformer. With conventional approaches, transformers connected to power plants are sized by assuming that all of the generators output their rated power without considering fluctuations which inherently characterize renewable power generation. It is then required to include this effect to have a more accurate sizing [4].

In the technical literature, rating of power transformers that connect wind power plants to the grid is investigated in [5], where a dynamic approach for exploiting the variability of power production is proposed. In [6] a cost optimization methodology is proposed which uses statistical scenarios based on distribution functions and Markov Chain models and accounts for maintenance and design requirements for transformer contingency in off-shore wind farms. The calculation of transformer limits according to the shape of the loading profile in presence of renewable generation is analyzed in [7], where an energy-limit concept for oil-immersed technology is proposed. Based on the loading factor, in [8] a method to size the power transformer using the thermal limit is proposed. The potential of dynamic rating in the planning stage of transformers for renewable sources is analyzed in [9] to utilize equipment more efficiently. In [10] a stochastic method is used to characterize daily loading profiles in the presence of fluctuations to evaluate their impact on the transformer lifetime. In [11] a rating method is proposed based on measured and estimated daily loading profiles. In [12] a data driven approach based on a Gaussian mixture mode is used to estimate the loading shape profile for the transformer rating at the purpose of long-term planning. Randomness of loading profiles in the presence of fluctuations are addressed in [13] through a Lognormal characterization of transformer lifetime for its power rating. The maximum utilization transformers through dynamic thermal rating is studied in [14] with reference to loading with random variations.

The analysis of the literature survey provides evidence of some critical issues that need to be addressed in the transformers’ rating, such as: (i) analysis of the loading conditions, (ii) thermal behavior of the transformer under specified loading conditions, (iii) estimation of the expected lifetime duration under specified loading conditions, and (iv) evaluation of the transformer size that better matches the desired lifetime and expected loading duration.

The analysis of the loading condition is clearly the main input of the design stage, where randomness, such as that regarding the fluctuation typical of wind farm applications, is faced with stochastic approaches, such as Markov Chain models [6], the Gaussian mixture [12], or the Wiener process [13]. In this regard, the proposed approach suggests and validates a novel rating method based on the statistical characterization of the loading condition through the Ornstein-Uhlenbeck process. This stochastic process appears particularly suitable for the description of the phenomenon under investigation for its ability in capturing the so-called ‘mean reversion’ of time series [15]. Based on the Ornstein-Uhlenbeck characterization of the loading current, the thermal description of the transformer is stochastically characterized as well. The solution of the thermal model, which is given by means of the Euler-Maruyama model, allows estimating lifetime reduction according to a Lognormal distribution. This is useful to derive an analytical method for the transformer sizing. Acceleration of lifetime reduction is avoided by selecting a transformer size larger than that obtained through the classical approaches based on average values of the loading current.

The paper is structured as follows: Section 2 describes the process used to stochastically characterize the wind power generation. Section 3 gives the details of the thermal model used for estimating the lifetime of a mineral-oil-immersed transformer and lifetime characterization in terms of stochastic processes. In Section 4 the analytical formulation of the transformer rating is calculated and the whole procedure is summarized. In Section 5 the results of the proposed rating method applied to the data of a wind farm are reported. Conclusions are drawn in Section 6.

2. Stochastic Characterization of the Wind Power

Many studies are available in the literature to properly characterize wind power randomness [16,17]. Among them, two approaches are typically used: static [18,19] and dynamic approaches [20,21]. The dynamic approach based on the stochastic process theory is the most suitable in order to achieve a more accurate characterization of the wind power behavior. As proposed in [20], the Ornstein-Uhlenbeck stochastic process appears particularly adequate to describe the random variations of this non-controllable energy source.

Indeed, the Ornstein-Uhlenbeck stochastic process is largely employed as a diffusion process to model several phenomena, as in the processing of financial high-frequency data. With respect to the power system, it is considered particularly suitable to describe the electrical load [22]. The rationale behind this choice relies in capturing the so-called ‘mean reversion’ of certain time series, which tend towards their mean values over the time. The Ornstein-Uhlenbeck process $P_t$, $t\ge 0$, is wholly described by the following stochastic differential equation:

$$d{P}_t=\tau \left(\mu -P_t\right)dt+\sigma d{W}_t$$

(1)

where µ is the long-term mean, τ is related to reversion speed, σ is linked to the instantaneous volatility, and $W_t$ is the Wiener process characterized by the well-known property that $W_t$ has independent increments with $W_t-W_s\approx N\left(0,t-s\right)$ (for $0\le t<s$), with $N\left(0,t-s\right)$ as the normal distribution with expected value 0 and variance t − s.

It is important to note that, if $0\le s_1<t_1\le s_2<t_2$, then $W_{t_1}-W_{s_1}$ and $W_{t_2}-W_{s_2}$ are independent random variables. Also, the following probabilistic properties can be derived:

-

the unconditional probability density function at a fixed time t is Gaussian:

$$f_{W_t}\left(x\right)=\frac{1}{\sqrt{2\pi t}}{e}^{\frac{-x^2}{2t}}$$

(2)

-

the expectation is zero:

$$E\left[W_t\right]=0$$

(3)

-

the variance at time t is:

$$Var\left(W_t\right)=E\left[W_t^2\right]-E^2\left[W_t\right]=E\left[W_t^2\right]=t$$

(4)

-

the covariance at time t is:

$$cov(W_s,W_t)=min(s,t).$$

(5)

The solution of the stochastic differential equation is given by:

$$P_t=P_0{e}^{-\tau t}+\mu \left(1-e^{-\tau t}\right)+\sigma {{\displaystyle \int }}_0^t{e}^{-\tau \left(t-s\right)}d{W}_t.$$

(6)

The Ornstein-Uhlenbeck process $P_t$ is stationary when $P_0\sim N\left(\mu ,\frac{{\sigma }^2}{2\tau }\right)$ and $P_0\perp W_t,t\ge 0$are assumed and has normally distributed increments and unconditional moments:

$$\begin{gathered} E\left[P_t\right]=\mu , \\ var\left[P_t\right]=\frac{{\sigma }^2}{2\tau },\mathrm{and} \\ cov\left[P_t,P_s\right]=\frac{{\sigma }^2}{2\tau }{e}^{-\tau \left|t-s\right|}. \end{gathered}$$

(7)

Non stationarity of the problem was attained for a given initial value $p_0$ with normally distributed increments and conditional moments:

$$\begin{gathered} E\left[P_t\lfloor P_0=p_0\right]=p_0{e}^{-\tau t}+\mu \left(1-e^{-\tau t}\right), \\ var\left[P_t\lfloor P_0=p_0\right]=\frac{{\sigma }^2}{2\tau }\left(1-e^{-2\tau t}\right), \\ cov\left[P_t\lfloor P_0=p_0\right]=e^{-\tau \left|t-s\right|}-e^{-\tau \left|t+s\right|},t\ne s. \end{gathered}$$

(8)

In order to investigate the problem of the process parameter estimation, the discrete version of the process is referred to here where, for the sake of simplicity, the time interval $T$ is assumed to be equal to 1. The process has $n$ observations $P_{t_i}$ at time instants $0=t_0<t_1<\dots t_n=T=1$, where $t_i-t_{i-1}=\frac{1}{n}$.

By referring to the discretized version of the process, the method of moments was applied for its simplicity. It is not difficult to demonstrate that the following relationships maintain as far as the estimate is concerned for the parameters μ, τ and σ² [22]:

$$\widehat{\mu }=\frac{1}{n+1}{\displaystyle \sum }_{i=0}^n{P}_{t_i}$$

(9)

$$\tau =n\log \frac{\frac{1}{n}{\sum }_{i=0}^n{\left(P_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)}^2}{\frac{1}{n-1}{\sum }_{i=1}^n\left(P_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)\left(P_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)}$$

(10)

$${\sigma }^2=2n\log \frac{1}{n}{\sum }_{i=0}^n{\left(P_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)}^2\log \frac{\frac{1}{n}{\sum }_{i=0}^n{\left(P_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)}^2}{\frac{1}{n-1}{\sum }_{i=1}^n\left(P_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)\left(P_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{P}_{t_i}\right)\right)}.$$

(11)

In practice, the actual observed process could be corrupted by a noise E_i; as a consequence of that, the following process should be examined:

$$Z_{t_i}=P_{t_i}+E_i.$$

(12)

Under the assumption E~N(0, ν²), we could obtain:

$$E\left[Z_t\right]=\mu ,var\left[Z_t\right]=\frac{{\sigma }^2}{2\tau }+{\nu }^2\mathrm{and}cov\left[P_i,P_j\right]=\frac{{\sigma }^2}{2\tau }{e}^{-\tau \left|t_i-t_{i-1}\right|}i\ne j.$$

(13)

In this case, the estimates for µ, τ, ${\sigma }^2$ and ν² can be given by the following relationships:

$$\widehat{\mu }=\frac{1}{n+1}{\displaystyle \sum }_{i=0}^n{Z}_{t_i}$$

(14)

$$\widehat{\tau }=n\log \frac{\frac{1}{n-1}{\sum }_{i=1}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}{\frac{1}{n-2}{\sum }_{i=2}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}$$

(15)

$${\widehat{\sigma }}^2=2n\frac{{\left(\frac{1}{n-1}{\sum }_{i=1}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\right)}^2}{\frac{1}{n-2}{\sum }_{i=2}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}\log \frac{\frac{1}{n-1}{\sum }_{i=1}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}{\frac{1}{n-2}{\sum }_{i=2}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}$$

(16)

$${\widehat{\nu }}^2=\frac{1}{n}{\sum }_{i=0}^n{\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}^2-\frac{{\left(\frac{1}{n-1}{\sum }_{i=1}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\right)}^2}{\frac{1}{n-2}{\sum }_{i=2}^n\left(Z_{t_i}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)\left(Z_{t_{i-1}}-\left(\frac{1}{n+1}{\sum }_{i=0}^n{Z}_{t_i}\right)\right)}.$$

(17)

3. Stochastic Thermal Model of a Mineral-Oil-Immersed Transformer

With lifetime evaluation being the core of the transformer’s design, the dynamic aspect of the thermal phenomena and their stochasticity due to the loading needed to be properly taken into account.

The thermal processes of mineral-oil-immersed transformers were modeled based on the top-oil and hot-spot temperatures. Top-oil temperature refers to the oil temperature at the top of the tank; hot-spot temperature refers to the temperature of the hottest point of the winding insulation. The measure of the top-oil temperature is simple and thus is suitable for on-line monitoring and control of the transformer. Hot-spot temperature is related to the loading capacity of the transformers and their lifetime degradation [23,24].

Regarding the top-oil temperature, the values it assumes depend on the heat generated by the Joule losses in the windings (i.e., load losses) and no-load losses, the ambient temperature, and the thermal characteristics of the oil. Regarding the hot-spot temperature, the values it assumes depends on the heat generated by the load losses, the temperature of the oil, and the thermal characteristics of both winding and oil [25,26,27,28,29]. In the case of a significant presence of non-linear loads or converter-based generators, the loads due to non-sinusoidal voltages can have a role in the thermal behavior of the transformer and consequently on its lifetime [30,31,32].

Estimation of both top-oil and hot-spot temperatures is typically made by thermal models which exploit the thermal-electrical analogy and heat transfer theory [1]. The thermal-electrical analogy theory is based on the use of simple circuits which allow analyzing the thermal behavior of the transformers through the similarity between thermal and electrical processes. Based on this approach, top-oil (${\vartheta }_{oil}$) and hot-spot (${\vartheta }_{hs}$) temperatures vary according to differential equations of resistive-capacitive lumped circuits forced by ideal sources related to heat generated in the transformer [1,33,34]. In Figure 1, the circuits related to the top-oil temperature (Figure 1a) and hot-spot temperature (Figure 1b) are reported.

In Figure 1a the total heat generated in the transformer ($q_{tot}$) refers to the heat generated by no-load losses ($q_{fe}$) and by the load losses ($q_l$), $C_{th-oil}$ is the equivalent thermal capacitance of the transformer oil, $R_{th-oil-air}$ is a nonlinear oil to air thermal resistance, and ${\vartheta }_{amb}$ is the known ambient temperature. In Figure 1b $q_{wdg}$ is the heat generated by winding losses, $C_{th-wdg}$ is the thermal capacitance of the winding, and $R_{th-hs-oil}$ is the nonlinear winding to oil thermal resistance. Non-linear thermal resistances in the circuits are considered to include the effect of the temperature on both oil thermal characteristics and Joule losses [24,33,34]. The non-linear thermal resistances allow obtaining more accurate temperature estimation due to the deep dependency of the oil characteristics—particularly the viscosity—on the temperature.

The differential equations governing the thermal circuits of Figure 1 are:

$$\frac{1+R{K}^2}{1+R}{\mu }_{pu}^{\nu }\mathsf{∆}{\vartheta }_{oil,rtd}={\mu }_{pu}^{\nu }{\tau }_{oil,rtd}\frac{d{\vartheta }_{oil}}{dt}+\frac{{\left({\vartheta }_{oil}-{\vartheta }_{amb}\right)}^{\nu +1}}{\mathsf{∆}{\vartheta }_{oil,rtd}^{\nu }}$$

(18)

$$K^2{P}_{cu,pu}\left({\vartheta }_{hs}\right){\mu }_{pu}^{\nu }\mathsf{∆}{\vartheta }_{hs,rtd}={\mu }_{pu}^{\nu }{\tau }_{wdg,rtd}\frac{d{\vartheta }_{hs}}{dt}+\frac{{\left({\vartheta }_{hs}-{\vartheta }_{oil}\right)}^{{\nu }^{\prime }+1}}{\mathsf{∆}{\vartheta }_{hs,rtd}^{{\nu }^{\prime }}}$$

(19)

with (18) being related to the top-oil model and (19) to the hot-spot model which depends on the top-oil temperature, with the following symbols meaning:

-

$R$ is the ratio of load losses at the rated current and no-load losses;

-

$K$ is the ratio of load current and the rated current, i.e., the load factor;

-

${\mu }_{pu}$ is the oil viscosity in per unit value;

-

${\vartheta }_{amb}$ is the ambient temperature;

-

${\vartheta }_{hs}$ is the hot-spot temperature;

-

${\vartheta }_{oil}$ is the top-oil temperature;

-

$\mathsf{∆}{\vartheta }_{oil,rtd}$ is the rated top-oil temperature rise over the ambient temperature;

-

$\mathsf{∆}{\vartheta }_{hs,rtd}$ is the rated hot-spot temperature rise over the top-oil temperature;

-

${\tau }_{oil,rtd}$ is the rated top-oil time constant;

-

${\tau }_{wdg,rtd}$ is the rated winding time constant;

-

$P_{cu,pu}$is the value of the load losses, in per unit, and $P_{cu,pu}\left({\vartheta }_{hs}\right)$ represents the dependence of these losses from the hot-spot temperature;

-

ν and ν′ are constant and depend on the type of transformer and the cooling modes; the exponentials of ν and ν$\prime$ determine the non-linearity of the model.

The solution of (18) and (19) can be evaluated through classical numerical methods, such as the Runge Kutta method, if a deterministic framework is considered. As will explained in the next section, the time variations of the top-oil and hot-spot temperatures have to be derived by properly treating the stochastic variations of the known load. For this purpose, the Euler-Maruyama method was employed.

The hot-spot temperature is the main factor responsible for accelerating the transformer aging, thus its estimation is used to evaluate the lifetime reduction of the winding insulation. Ambient temperature and moisture also have influence, the latter having a significant effect since the water generated by the process of cellulose decomposition significantly accelerates its aging rate [35,36,37].

The insulation deterioration in the time period $T$ can be evaluated on the basis of experimental evidence according to [38]:

$${\eta }_M=\frac{1}{T}\underset{0}{\overset{T}{\int }}{F}_{AA}dt$$

(20)

with ${\eta }_M$ being the loss of life of the transformer over the considered time period $T$ and $F_{AA}$ being the aging acceleration factor. Regarding the duration of the time period to be considered for designing purposes, it can be assumed that $T=24\mathrm{h}$, since the working cycle is typically repeated according to daily working cycles [39]. Regarding the aging acceleration factor $F_{AA}$, it is governed by the well-known Arrhenius reaction rate theory, which is given by:

$$F_{AA}=e^{\frac{\mathrm{15,000}}{110+273}-\frac{\mathrm{15,000}}{{\vartheta }_{hs}+273}}$$

(21)

where the hot-spot temperature ${\vartheta }_{hs}$ is assumed time-varying and 110 °C is assumed as the reference temperature, i.e., the acceleration factor of degradation is evaluated based on the deviation of ${\vartheta }_{hs}$ from 110 °C. In this way, the lifetime decrease ${\eta }_M$ was estimated through (20) in per unit measurements.

3.1. Integration of the Stochastic Differential Model

By neglecting the voltage drop across the leakage reactance, the load factor $K$ in (18) and (19) can be considered equal to $\left(\frac{P_t}{S_n}\right)$, with $S_n$ being the base power. Due to its dependence on the process $P_t$, $K$ is properly described as random process $K_t$, as well as the ratio of load losses$R$, which is also described as random process $R_t$. For characterizing the stochasticity of the variables ${\vartheta }_{oil}$ and ${\vartheta }_{hs}$, the following model had to be integrated:

$$d{P}_t=\tau \left(\mu -P_t\right)dt+\sigma d{W}_t$$

(22)

$$d{\vartheta }_{{oil}_t}=\frac{{\left({\vartheta }_{{oil}_t}-{\vartheta }_{amb}\right)}^{\nu +1}}{{\mu }_{pu}^{\nu }{\tau }_{oil,rtd}\mathsf{∆}{\vartheta }_{oil,rtd}^{\nu }}dt-\frac{1+R_t{K}_t^2}{{\tau }_{oil,rtd}\left(1+R_t\right)}\mathsf{∆}{\vartheta }_{oil,rtd}dt$$

(23)

$$d{\vartheta }_{{hs}_t}=\frac{{\left({\vartheta }_{{hs}_t}-{\vartheta }_{{oil}_t}\right)}^{{}^{\nu \prime }+1}}{\mathsf{∆}{\vartheta }_{hs,rtd}^{{}^{\nu \prime }}{\mu }_{pu}^{\nu }{\tau }_{wdg,rtd}}-\frac{K_t^2{P}_{{cu,pu}_t}\mathsf{∆}{\vartheta }_{hs,rtd}dt}{{\mu }_{pu}^{\nu }{\tau }_{wdg,rtd}}.$$

(24)

The integration of this model was numerically performed by using the Euler-Maruyama method, which is opportune to briefly recall [40]. Let us consider a scalar, autonomous stochastic differential equation written in the standard form:

$$dX\left(t\right)=g\left(X\left(t\right)\right)dt+h\left(X\left(t\right)\right)d{W}_{t}X\left(0\right)=X_{0}0\le t\le T.$$

(25)

The discretized version of (25) is

$$X_i=X_{i-1}+g\left(X_{i-1}\right)\mathsf{∆}t+h\left(X_{i-1}\right)\left(W_t{}_i-W_t{}_{i-1}\right)i=1,\dots \dots .,n$$

(26)

where $\mathsf{∆}t=\frac{T}{n}$.

In (26) the only complication lies in the generation of the increments $\left(W_t{}_i-W_t{}_{i-1}\right)$. By keeping in mind the properties of the standard Wiener process, also named standard Brownian motion, it is intuitive to assume $\left(W_t{}_i-W_t{}_{i-1}\right)$ is an independent random variable distributed according to the law $\sqrt{\Delta t}N\left(0,1\right)$. The Euler-Maruyama method is substantially summarized by Formula (26), where only the additional burden to generate random vector from the $N\left(0,1\right)$ distribution is required.

It is good to make some considerations regarding convergence of this method, that is, to ascertain whether the solution tends more and more to true when $\mathsf{∆}t$ tends to zero. By denoting with $X_i^{*}$ the true solution at the ith-step, which is a random variable as $X_i$, the concept of strong convergence has to be introduced. A method exhibits strong order of convergence equal to α if a constant λ exists such that:

$$E\left[\left|X_i-X_i^{*}\right|\right]\le \lambda \Delta t^{\alpha }i=1,\dots \dots .,n.$$

(27)

It can be demonstrated that if the functions $g$and $h$ satisfy some conditions, the Euler-Maruyama method has strong order of convergence $\alpha =1/2$. The concept of strong convergence is stated in terms of expected value, but useful inequality can also be derived to characterize the error for individual simulations by exploiting the Markov inequality. As is well known, this cornerstone inequality in the theory of probability states that if a random variable $X$ has finite expected value, $\forall a>0$:

$$P\left(\left|X\right|>a\right)\le \frac{E\left|X\right|}{a},i=1,\dots \dots .,n.$$

(28)

By applying this inequality to (27), and choosing $a=∆t^{\mathsf{\delta }}$, with $\mathsf{\delta }<\alpha$, since for Euler-Maruyam method $\alpha =1/2$, it can be written as:

$$P\left(\left|X_i-X_i^{*}\right|>\Delta t^{\delta }\right)\le \frac{E\left[\left|X_i-X_i^{*}\right|\right]}{\Delta t^{\delta }}\le {\lambda \Delta t}^{\left(\frac{1}{2}-\delta \right)}.\hspace{1em}\hspace{1em}i=1,\dots \dots .,n$$

(29)

It can be more useful to write the relationship as:

$$P\left(\left|X_i-X_i^{*}\right|<\Delta t^{\delta }\right)\ge {1-\lambda \Delta t}^{\left(\frac{1}{2}-\delta \right)}\hspace{1em}\hspace{1em}i=1,\dots \dots .,n$$

(30)

which ensures that for any point over the interval $\left[0,T\right]$, the error will be small with the probability tending to $1$.

3.2. Lifetime Characterization

The oil transformer temperature as well as the hot-spot temperature have the character of stationary random variables. As will be numerically demonstrated in the numerical application, the mean value and the standard deviations can be considered as constant over the whole interval $T$. More specifically the hot spot temperature can be optimally approximated by a stationary Gaussian process.

The integrand factor in Equation (20) can be rewritten as:

$$F_{AA}\left(t\right)=e^{\frac{\mathrm{15,000}}{110+273}-\frac{\mathrm{15,000}}{{\vartheta }_{hs}\left(t\right)+273}}=e^{k^{\prime }\Delta {\vartheta }^{\prime }\left(t\right)}$$

(31)

with

$$k^{\prime }=\frac{\mathrm{15,000}}{{\left(110+273\right)}^2}$$

(32)

$$\Delta {\vartheta }^{\prime }\left(t\right)={\vartheta }_{hs}\left(t\right)-110.$$

(33)

$F_{AA}$ is a Lognormal process and its parameters can be expressed explicitly as a function of the moments of the stochastic process, $∆{\vartheta }^{\prime }$. The Lognormal statistical feature of $F_{AA}$ and its stationarity were verified in the numerical application reported in Section 5 by examining the behavior of both the mean value and the standard deviation over the time.

From the earlier definition of the reduction of the transformer’s lifetime, the following relationship can be stated:

$${\eta }_M=\frac{1}{T}{\int }_0^T{F}_{AA}dt=\frac{1}{T}{\int }_0^T{e}^{k^{\prime }\mathsf{\Delta }{\vartheta }^{\prime }}dt.$$

(34)

As a consequence of the stationarity of the $F_{AA}$ process, the reduction in the transformer’s lifetime ${\eta }_M$can be optimally approximated by a Lognormal random variable whose mean value is equal to

$$E\left[{\eta }_M\right]=e^{k^{\prime }{\mu }_{\mathsf{\Delta }{\vartheta }^{\prime }}+\frac{k^{\prime 2}{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}}.$$

(35)

When $E\left[{\eta }_M\right]=1$ (that is $k^{\prime }{\mu }_{\mathsf{\Delta }{\vartheta }^{\prime }}+\frac{k^{\prime 2}{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}=0$) there is no ageing acceleration with respect to that arising with the rated values of current, $I_{rtd}$. The fluctuation of the loading current around its mean value $E\left[I_l\right]$ (that is the value typically used to size transformers $I_{rtd}=E\left[I_l\right]$), however, implies an ageing acceleration (i.e., lifetime reduction) higher than that occurring for $I_{rtd}$ and whose value is related to $k^{\prime }{\mu }_{\Delta {\vartheta }^{\prime }}+\frac{k^{\prime 2}{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}\ne 0$ corresponding to the variation of temperature due to the current’s fluctuation.

4. The Proposed Transformers Rating Procedure

The value of $I_{rtd}$ (i.e., of transformer’s size) can be identified by assuming that $E\left[I_l\right]$ produces a constant temperature equal to $\mathsf{∆}{\vartheta }_{hs,rtd}-{\mu }_{\mathsf{∆}{\vartheta }^{\prime }}-\frac{k^{\prime }{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}$ and the fluctuating current the contribution ${\mu }_{\mathsf{∆}{\vartheta }^{\prime }}+\frac{k^{\prime }{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}.$ By keeping this in mind, the transformer rated current can be estimated as follows:

$$I_{rtd}=E\left[I_l\right]\gamma$$

(36)

with

$$\gamma =\sqrt{\frac{\mathsf{∆}{\vartheta }_{hs,rtd}}{\mathsf{∆}{\vartheta }_{hs,rtd}-{\mu }_{\mathsf{∆}{\vartheta }^{\text{’}}}-\frac{k^{\prime }{\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}^2}{2}}}.$$

(37)

From (37), it is possible to straightforwardly obtain the rated value of the current $I_{rtd}$. Indeed, the design procedure has to start by initially choosing a transformer whose rated power is as near as possible to the mean value of the loading power. By making reference to its characteristic parameters, one can perform an a priori estimation of the temperature distribution. This allows directly determining that the oversizing factor $\gamma$ for ensuring the mean value of lifetime is equal to that expected.

The procedure for choosing the rated power of the transformer can be eventually summarized by the flow chart reported in Figure 2.

The proposed procedure—shown in Figure 2—can be summarized as follows:

• the actual power profile of the wind farm is recorded based on data available from measurements;
• the actual power profile is statistically analysed as an Ornstein-Uhlenbeck process and a number of synthetic profiles are generated;
• based on the loading current corresponding to the synthetic power profiles, the thermal model of the transformer is applied and solved in statistical terms through the Euler-Maruyama method, thus obtaining the statistical characterization of the hot-spot temperature;
• on the basis of the distribution of the hot-spot temperature, the statistical characterization of the transformer’s lifetime is derived;
• the rating of the transformer able to guarantee the desired lifetime is calculated through Equation (37).

5. Numerical Application

In this section, the results of the application of the proposed rating method to the data of a 600 MW wind farm are reported [41]. Based on the measured data (step 1 of the flow chart in Figure 2), the actual power produced by the wind farm is plotted in Figure 3 (blue thicker line). This dataset refers to the time horizon of one day which is divided into 96 time-intervals, each 15 min long.

This daily profile is characterized as an Ornstein-Uhlenbeck process (step 2 of the of the flow chart in Figure 2), whose parameters, estimated by the method of moments, are $\tau =25.56$, $\mu =506.40$, and $\sigma =530$. Some realizations of the process (i.e., synthetic profiles), by employing these parameters, are also reported in Figure 3 (gray lines).

The rating of the transformer connecting the wind farm to the network was performed according to the method proposed in this paper. The parameters used in the thermal model were derived from a typical 605 MVA OFAF transformer [33] and are reported in

Table 1. Parameters for the thermal model of the transformer.

Parameter	System	Value
Rated heat generated by load losses	$q_{cu}$	930 kW
Rated heat generated by no-load losses	$q_{fe}$	285 kW
DC loss per unit value	$P_{cu,pu}$	0.9
Eddy losses per unit value	$P_{cu,eddy,pu}$	1.1
Ambient temperature	${\vartheta }_{amb}$	30 °C
Rated top-oil temperature rise over ambient temperature	$\mathsf{∆}{\vartheta }_{oil,rtd}$	35 °C
Rated hot-spot temperature rise over top-oil temperature	$\mathsf{∆}{\vartheta }_{hs,rtd}$	45 °C
Heat exchange area	$A$	2580 m2
Empirical value constants	$C$	0.59
	ν	0.25
Characteristic dimension	$L$	10 m
Oil thermal conductivity	$k$	0.128 W/mK
Oil mass	$m_{oil}$	79,746 kg
Winding time constant	${\tau }_{wdg}$	12 min

. The ambient temperature was assumed to be constant (30 °C), since the assumption was made that a confined ambient environment equipped with a cooling system would almost guarantee a constant temperature. Regarding the oil thermal conductivity, its variation with the temperature was negligible compared to the variations of the oil viscosity [32]. Thus, the approximation of constant thermal conductivity corresponding to 65 °C was assumed. Both C and ν are empirical constants whose values correspond to the case of laminar oil circulation [32].

In order to perform the proposed method for the choice of the transformer rating, 10,000 daily profiles of wind power were obtained through the statistical characterization of the Ornstein-Uhlenbeck approach. These profiles were used to derive the distribution of hot-spot temperature at each time interval of the day. With reference to the rated currents corresponding to the mean values of this data set, in Figure 4 and Figure 5 the mean value and standard deviation of this temperature at each interval are reported. The mean values reported in Figure 4 reveal that they vary within a narrow interval around the rated hot-spot temperature, i.e., 110 °C, which corresponds to a standard lifetime duration assumed for the transformer. The profile of the standard deviation at each interval reported in Figure 5 clearly converges to a specific value, thus showing that the considered process can be assumed to be stationary.

Regarding the hot-spot temperature, it can be stated that its distribution at each interval is Gaussian. This is verified by the results obtained by the adoption of the thermal model of Section 3 (step 3 of the of the flow chart in Figure 2). Particularly in Figure 6, the histogram of the hot-spot temperature corresponding to a generic interval of the day (the interval corresponding to the 12th hour) is reported together with its fitting Gaussian probability density function (pdf) which shows a good approximation with the simulated data. Mean value and standard deviation of this distribution are ${\mu }_{{\vartheta }_{hs}}=110.92$ °C and ${\sigma }_{{\vartheta }_{hs}}=5.91$ °C.

The accurate Gaussian approximation shown in Figure 6 is coherent with the assumption of the Ornstein-Uhlenbeck model. The analysis of the other intervals of the day, not reported here for the sake of brevity, shows that this Gaussian approximation is also valid.

The distribution of the aging acceleration factor at the same interval of Figure 6 is reported in Figure 7.

As expected, this distribution is well fitted by a Lognormal distribution, whose pdf is also reported in the figure. The analysis of the other intervals of the day, here not reported for the sake of conciseness, shows that this Lognormal approximation is also valid. This approximation is coherent with (31) and with the Gaussian approximation of the hot-spot temperature at each interval of the day (Figure 6).

As is well known, the Lognormal distribution can be represented in its general form as:

$$f\left(x\right)=\frac{e^{\frac{-{\left(\ln x-\mu \right)}^2}{2{\sigma }^2}}}{x\sqrt{2\pi }\sigma }$$

(38)

where $\mu$ and $\sigma$ are the mean and standard deviation of the associated Gaussian distribution ($\mu =0.10$ and $\sigma =0.65$in our case).

Based on the Lognormal process of the aging acceleration factor, the distribution of lifetime reduction in per unit values was derived and is shown in Figure 8 (step 4 of the of the flow chart in Figure 2).

Coherently with the theoretical hypothesis made in Section 3, it was found that the distribution of the lifetime reduction values is well approximated by a Lognormal pdf, as clearly appears in Figure 8. This can be explained by considering that the lifetime is evaluated by means of (34) and the aging acceleration factor (31) is well approximated by a Lognormal distribution at each interval (Figure 7). In this case, even if the sum of Lognormal distributions is not generally lognormally distributed, the stationary property of the process implies that the Lognormal distribution can still be used to represent the integral in (34). It is worth noting that this result has general validity, since it derives from the statistical characterization of the wind power according to the Ornstein-Uhlenbeck stochastic process and a consolidated thermal model of transformers. In order to apply the proposed procedure at the design stage of the transformer with the purpose of choosing the rating of the transformer, the proposed approach was repeatedly applied by increasing values of the rated transformer current, as recalled in the flow chart of Figure 2 (step 5 of the flow chart in Figure 2). Figure 9 shows the lifetime reduction corresponding to the ratio $\left(\frac{I_{rtd}}{E\left[I_l\right]}\right)$ of the rated and mean values of the loading current.

Figure 9 clearly shows that the reduction of lifetime—compared to the duration corresponding to the standard hot-spot temperature (110 °C) decreases when the rated current increases. Particularly, with reference to the proposed case study, the lifetime reduction is equal to one when the rated current is 1.032 times greater than the mean value of the loading current. For lower values of the rating current the reduction of lifetime is greater than one, reaching a value higher than 1.25 in case of rated current equal to the mean value of the loading current. This means that by applying conventional rating methods, i.e., by neglecting the effect of fluctuations of the generated power, the expected lifetime of transformer would be reduced by 25%. Alternatively, to avoid this reduction, a de-rating of the transformer would be needed. Figure 9 also shows that in the case of higher values of rated current (i.e., abscissa greater than 1.032) a slow reduction of lifetime appears (i.e., its values are lower than one), reaching the value of 0.6 when the rated/loading current ratio is 1.1.

In order to analytically derive the increased rated current identified in Figure 9, Equation (36) was also applied with ${\mu }_{\mathsf{∆}{\vartheta }^{\prime }}=0.93°\mathrm{C}$ and ${\sigma }_{\mathsf{∆}{\vartheta }^{\prime }}=5.91°\mathrm{C}$. Coherently with the result of Figure 9, the rated current corresponding to the rated hot-spot temperature was equal to 1.032. This result once again shows the validity of the hypothesis that lifetime reduction is lognormally distributed, as derived in Figure 8.

6. Conclusions

In this paper, the optimal rating of transformers connecting wind power plants has been dealt with by taking into account the probabilistic analysis of lifetime degradation. The Ornstein–Uhlenbeck stochastic approach was used to generate synthetic wind turbines power profiles, as well as synthetic thermal transformer profiles which allowed formulating the aging of transformers. In the numerical application, the lifetime reduction was demonstrated to be distributed as a Lognormal process allowing to easily derive the rated current corresponding to an aging acceleration factor equal to that obtained in correspondence of the rated hot-spot temperature, i.e., the standard lifetime duration assumed for the transformer. The results of the simulation showed that the proposed method is able to determine the increment of rated current required to get a standard lifetime duration without incurring a higher aging acceleration. This is particularly useful in the case of wind farm applications whose high wind power fluctuations do not modify the mean power value which leads to incorrect transformer sizing. The use of Ornstein–Uhlenbeck stochastic approach is able to capture the effect of these fluctuations in terms of lifetime reduction. Future application will be devoted to the combination of the proposed method with other sizing issues, which typically affect the aging of transformers of wind farms, including through the harmonic content, due to the use of wind power converters.