Bayesian Inference for Thermal Model of Synchronous Generator—Part I: Parameter Estimation

Due to the increasing injection of intermittent power sources (solar+wind) into a common grid, dispatchable sources such as hydro power should be able to help reduce the variability in load and the variability in generation caused by the intermittent sources. A hydro generator should be able to operate short-term beyond its thermal capability limit. This requires the monitoring of internal temperatures in the hydro generator. In this paper, a thermal model of an air-cooled synchronous generator is presented, emphasizing the various aspects of parameter estimation and identifiability using Bayesian inference. Inferences are drawn from the posterior distributions of the parameters and initial conditions, dispersion (spreading) of particles and sampling efficiency, practical parameter identifiability, and model mismatch with experiments. Results show extremely narrow parameter distributions. It is early to generalize about the posterior distribution of air-related and metal-related parameters of the air-cooled synchronous generator based on the single experimental data presented here.


I. INTRODUCTION
Electricity generation from intermittent sources such as solar power, wind power, tidal power, etc., is rapidly increasing in modern electric power system networks. The intermittency in these sources causes the power system networks to operate in different operating conditions. Dispatchable sources such as hydro power can be used for removing the variability in the system's power production caused by the intermittent sources [1], [2], [3], [4]. Thus, in a modern power system, the hydro generators play a significant role in the flexible operation of the intermittent grid. A concept of flexible hydro power is coined in [5] for modern intermittent power system networks. This adheres to a new research requirement in the case of a synchronous hydro generator operating in tandem with the intermittent sources. The performance of the synchronous generator depends on its capability diagram [6]. The capability diagram provides information about the oper-The associate editor coordinating the review of this manuscript and approving it for publication was Bo Shen . ating regimes of the synchronous generator in case of the various operational limits, viz., armature current limit, field current limit, and under-excitation [7]. In [6], an instance of exploiting more active power from the hydro generator is studied by controlling the internal temperature of the machine. The control of the armature current limit and the field current limit will result in a decrease in resistance of the armature and the field winding due to temperature monitoring of the rotor copper, stator copper, and stator iron. Furthermore, because of an increase in the active current through the synchronous generator, more active power can be exploited. The temperature is controlled by cooled air circulation through the generator's internal surfaces. The cooled air is supplied through a heat exchanger, in a closed loop. Previous work includes a brief review of thermal analysis of electrical machines [8]. Lumped-parameter thermal network (LPTN) models of the thermal machines are provided in [9], [10], and [11]. Finite element analysis (FEM), and computational fluid dynamics (CFD) models were studied in [12] and [13]. A totally enclosed water-cooled thermal model 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/ of synchronous machines for an electric vehicle has been proposed in [14]. More recently, a totally enclosed thermal model of an air-cooled hydro generator has been developed in [15] using a closed-loop heat exchanger model for cooling heated air from the outlet of the generator. The thermal model of the air-cooled generator is further extended in [16] with the inclusion of temperature-dependent electrical resistances, temperature-dependent specific heat capacities of the metals, and fluids (air+water) inside the air-cooled hydro generator. The model of the air-cooled hydro generator is represented by a computationally cheap online solution of the non-linear model of the heat exchanger in [17], where a hybrid model (mechanistic+data-driven) is proposed using linear and nonlinear regression.
In Section II, materials and methods are outlined. Section III provides the mathematical governing equations for the air-cooled hydro generator. In Section IV, results from parameter estimation and parameter identifiability using Bayesian inference are discussed. Conclusions are drawn in Section V with future work suggested in Section VI. Figure 1 shows the working principle of a thermal model of an air-cooled synchronous hydro generator. The cold air out of the heat exchanger is blown by a fan into the rotor/stator air gap. The air is heated by heat flow from the rotor, air gap windage, and bearing friction. Furthermore, the air is forced into the iron cores and then gets heated by the heat flow from the iron cores. The heated air is now collected at the stator's outlet and passed through a counter current heat exchanger. The heated air is cooled using continuous cold water circulation in the heat exchanger and then fed again into the air gap as a closed loop process. The heat exchanger is provided with cold water, with mass flow ratė m w , at temperature T c w . The air mass flow rate isṁ a with temperature T h a at stator outlet and heat exchanger entry. The rotor copper heat source,Q σ r , is due to eddy currents caused by I f . Similarly, the stator copper heat source,Q σ s , is due to stator terminal current I t .Q σ Fe is stator iron heat source, andQ σ f is heat generated due to friction in the stator/rotor air gap. The thermal operation of the air-cooled synchronous generator is mainly influenced byṁ w ,ṁ a , T c w ,Q σ Fe ,Q σ f , I t , and I f . It is of interest to consider evolution of temperature in the rotor, stator, and iron core indicated by T r , T s and T Fe , respectively. Monitoring of these temperatures allows for optimal exploitation of active power production by enhancing the capability diagram to a new regime of operation [6]. Figure 2 shows the Bayesian framework for inference about parameters of a dynamical system. In the figure, x, u, z, θ, and y are the states, inputs, algebraic variables, parameters, and outputs, respectively. In the figure, p (θ) is the prior probability distribution of θ, p (y | θ) is the likelihood function, and p (θ | y) is the estimated posterior distribution of θ. Section IV provides detailed explanation about  priors and likelihood. In Fig. 2, we also see that the posterior distribution of parameters allows for various inferences about the parameters of a dynamical system. The inferences include (i) finding statistical moments from the posterior distribution of the parameters, (ii) finding the relative dispersion of the parameter space and the relative sampling efficiency between the parameters, (iii) model mismatch with experiment and (iv) inferences related to the posterior distribution of the initial conditions while working with the model and the experimental data offline.

B. EXPERIMENTAL DATA
A heat-run test, of the synchronous machine, was performed for 600 min [15] with sampling rate = [1] min. Only data from t = 16 min to t = 600 min, i.e., 584 data points, will be used for model fit. For each minute, for a supplied field current, starting from a cold-start, measurements for different quantities are recorded. The cold-run last up to 53 min, where the terminal voltage is built-up due to residual flux in rotor windings. The measurements are available for both electrical quantities and temperatures related to the aircooled synchronous generator. After the cold-run, the field current is increased which increases the temperature of the stator copper and the stator iron. The measured quantities are summarized in Table 1, and the experimental data are plotted in Fig. 3. The expression for terminal current I t as shown in Table 1, indicates that it is not measured using a sensor, however calculated from a mathematical expression relating power and voltage.

III. MATHEMATICAL MODEL
The mathematical equations governing generator metal temperatures taken from [16] and [18] are  Experimental data for generator model from a 600 min heat-run test taken from [15].
Similarly, the dynamical equations for air inside the generator are and the heat exchanger is model as St for i ∈ {w, a, } are dimensionless Stanton numbers relating heat transfer coefficient, density, heat capacity, and velocity.
Equations 1-6 can be written in Differential Algebraic Equations (DAEs) form as . The parameters and operating conditions are given in Table 2.
Out of the three states, T s and T Fe are measured, while it is of interest to estimate the temperature of rotating rotor copper T r . Similarly, out of three algebraic variables, T c a and T h a are measured, and it is also of interest to estimate air gap temperature T δ a . The measured inputs, states and algebraic variables are shown in Fig. 3.

IV. PARAMETER ESTIMATION A. PROBLEM FORMULATION
It is of interest to estimate the thermal parameters and initial conditions of the air-cooled hydro generator using Bayesian inference.
The expected value of parameterθ is calculated aŝ where p (θ | y) is the ''posterior probability distribution of parameter θ for the given data y''. p (θ | y) is expressed in terms of likelihood p (y | θ) and prior p (θ), where p (y) is independent of θ, and is used as a normalization constant for p (θ | y). p (y) is also known as the evidence or the marginal likelihood. The prior p (θ) is our prior beliefs about the probability density function for the parameter θ without seeing the data. Similarly, the likelihood p (y | θ) is a model representing the distribution of the data given a fixed parameter θ and calculated as The evidence p (y) used for normalization is calculated as the joint probability distribution of p (θ | y) and p (y) The analytical solution to Eq. (9) is only available for simple cases, and in real-life Bayesian inference, numerical solutions are used. The Bayesian parameter estimation method is applied to the system given by Eqs. (1-6).

B. FORMULATION USING TURING.JL
Turing.jl is a Julia package for probabilistic programming [19]. Turing.jl is also composable with DifferentialEquations.jl [20], a Julia package for differential equations, facilitating Bayesian inference in the parameter of the system VOLUME 10, 2022  represented by differential equations. The problem is formulated as shown in Fig. 4. For the differential equations represented by Eqs. (1-6) we want to estimate posterior distributions of: Table 3 shows the priors and the data for the estimation of the selected parameters. In the table, V, usually chosen as a inverse gamma function mostly for measurement sequence in Bayesian inference [21], is the prior to the variance of the measurement noise where we assume V ∼ −1 (2, 3). We have assumed that all measurements have the same noise variance since all temperatures are of similar size. The measurement vector is y = T s T Fe T c a T h a . A loss function is formulated using Turing.jl. The priors for the initial rotor copper temperature are truncated normal distributions. For parameters θ such as initial descriptor (differential and algebraic variables) as well as model constants, it is common to assume a normal distribution, e.g., θ ∼ N (µ, σ ). Because we normally want to limit the distribution to lie within a range θ ∈ [θ min , θ max ], e.g., to avoid negative values, it is quite common to use a truncated normal distribution for the prior of θ, where T represents truncated normal distribution for θ. The prior to the initial rotor copper temperature is chosen as T r (t = 0) ∼ T (N (30, 3) , 25, 35) where the mean µ T r (t=0) = 30 • C and is taken from Table 2 and the standard deviation σ T r (t=0) = 3 is the initial deviation that is assumed. Similarly, the prior for T r (t = 0) is truncated between 25 • C and 35 • C. The priors of T s (t = 0) and T Fe (t = 0) are also set accordingly with variance = 3 and the mean value taken from Table 2 within some relevant values of θ min and θ max for the parameters. Priors of other parameters to be found are also set accordingly from Table 2. It is important to note that the posterior distributions are approximated based on the numerical solution of Eq. (9) using different sampling methods. It is out of the scope of this paper to detail sampling methods. Some of the available sampling methods are listed in Table 3 and usage of these sampling algorithms can be found in [19]. We have chosen the NUTS sampler with the number of particles in the sampling as N s = 1000 to estimate the parameters.

C. ESTIMATED PARAMETERS
The posterior distributions of the estimated heat transfer parameters are shown in Fig. 5. Figure 5 (a) shows the distribution of the variance V of the measurement noises.    Table 4. The expected value of these distributions are used as the estimated parameters. Similarly, Fig. 6 shows the estimated initial conditions for the metal temperatures of the air-cooled generators.

D. ANALYSIS OF ESTIMATED PARAMETERS
The estimated parameters are analyzed based on mean and standard deviation as well as naive standard error (Naive SE) [22], and effective sample size (ESS) [23] as shown in Table 4.

1) NAIVE STANDARD ERROR (NAIVE SE)
Naive SE is a term defined for inferential statistics similarly to the mean and the standard deviation defined for descriptive statistics. Naive SE is computed as in [22]. Naive SE provides a measure of the potential error in the estimate while parameter inference is done through Bayes' theorem. Naive SE  calculates the width of sample means around the population mean. The lower the value of naive SE, the lower is the dispersion. Naive SE can also be used to find the upper and lower limit for the 95% confidence interval of the parameter given byθ ±σ θ whereθ is the mean value of the parameter andσ θ is the naive SE.
From Table 4 we see that Naive SE measuring the dispersion of sample means around the population mean, in the case of initial conditions lies between 4.74 · 10 −14 to 6.48 · 10 −14 .

VOLUME 10, 2022
The ratios of Naive SE for initial values isσ T r :σ T s :σ T Fe = 4.74 : 6.48 : 6.12. Since T r (t = 0) has lower Naive SE, the posterior distributions of T s (t = 0) and T Fe (t = 0) are wider than the posterior distribution of T r (t = 0). This shows that the posterior distribution of the temperature of the initial states related to the rotating copper inside the hydro generator is less wider than that of stationary copper and iron.
From Table 4 we see that the naive SE in the case of heat transfer from rotor copper to the air-gap UA r2δ is smaller as compared to the heat transfer from the stator copper to stator iron UA s2Fe and the heat transfer from iron to air UA Fe2a . Thus, the posterior distribution of the heat transfer parameter related to the air is more narrower than the heat transfer parameter related to the metal inside the hydro generator.
From the values of the naive SE from Table 4 in the case of heat sources parameters of the hydro generator, the dispersion in the posterior distribution of the iron heat source parameteṙ Q σ Fe is higher than that of the friction heat source parameteṙ Q σ f indicating that the heat source parameter related to air is less wider than the metals.

2) EFFECTIVE SAMPLE SIZE (ESS)
ESS describes the correlation between observations in the sample [23]. The calculation of ESS in Turing.jl is performed as in [24]. The higher the value of ESS, the higher is the correlation between the observations in the sample. ESS also helps to determine the relative sampling efficiency of the estimates based on the correlation between the observations in the sample. The relative sampling efficiency of the parameters from the estimated ESS can be calculated aŝ The smaller the value of ESS, the higher is the sampling efficiency. From Table 4, the relative sampling efficiencŷ η among the initial conditions T r (t = 0), T s (t = 0), and T Fe (t = 0), the relative sampling efficiency of stator iron temperature T Fe (t = 0) is higher as compared to the relative sampling efficiency of T r (t = 0) and T s (t = 0). In addition comparing all the temperatures, the rotor copper temperature T r (t = 0) has the lowest sampling efficiency. This shows that the temperature related to the rotating part of the hydro generator has lower sampling efficiency than the stationary part of the hydro generators. The sampling efficiency of the temperatures can be calculated using Eq. (10). The sampling efficiency for rotor copper temperature is given asη T r = 1 −  Table 4 using the estimated ESS for the heat transfer parameters, the relative sampling efficiency of heat transfer from rotating rotor copper to air-gap UA r2δ is lower as compared to heat transfer between stationary copper to stationary iron UA s2Fe or heat transfer from stationary iron to air UA Fe2a . Similarly from Table 4, we can see that both heat source parameters have the same sampling efficiency.

E. PARAMETER IDENTIFIABILITY
Parameter identifiability tells whether a parameter can be computed uniquely from the given model structure and observations. For complex systems, the number observed quantities is much smaller than the number of states + algebraic variables. It is therefore of interest to estimate the distribution of parameters that can explain the experimental data well. In inferential statistics, the joint posterior distribution of parameters found in Section IV using the Bayesian method can be used for parameter identifiability analysis. In frequentist statistics, profile likelihood projections are used for parameter identifiability [25], [26], [27]. Parameter non-identifiability occurs due to (i) indistinguishability of parameters in the model structure, and (ii) insufficiency in the experimental data. Identifiability analysis considering model structure is termed structural identifiability and identifiability analysis considering the amount and quality of experimental data is termed practical identifiability. Structural identifiability is out of the scope of the paper and our focus is on practical identifiability analysis. Figure 7 shows three cases of the joint posterior probability distribution of parameters θ 1 and θ 2 where Fig. 7 (a) illustrates that both parameters are structurally non-identifiable since the parameters do not converge to a point. The white lines show the posterior high density interval (HDI) within which an unobserved parameter value falls with a particular probability. The white dashed line indicates that the parameters diverge to infinity. The non-identifiability in parameters can only be resolved after the model structure is distinguishable with parameters. In Fig. 7 (b), the parameters are partially identifiable only at the lower density interval. The partially identifiable parameters are denoted as practically non-identifiable parameters. Identifiability of practically non-identifiable parameters can be improved by increasing the amount and the quality of the experimental data. Finally, Fig. 7 (c) shows that parameters converge to a point and are identifiable. Figure 8 shows the posterior joint probability distribution or the marginal kernel density estimate of heat transfer parameter UA r2δ with other heat transfer and heat source parameters. The central region in the marginal posterior distribution shows the HDI for the parameter space with a higher confidence region for the estimated parameters. Since the joint density plot of other heat transfer and heat source parameters bounded within a region, all the heat transfer and heat source parameters are identifiable. The joint points in the central region bounded with distorted ellipses in the figure show the The identifiability of all the parameters and initial conditions can also be inferred directly from the posterior distribution of the parameters shown in Figs. 5 and 6, respectively.

F. MODEL FITTING
It is of interest to see how well the mathematical model fits with the experimental data. The estimated initial conditions and parameters are used to compare the fitted model with the experimental data. Figure 9 shows the simulation versus experiment using the estimated parameters for the simulation. The model represents the experimental data well. In the figure, experimental data are well-matched with the simulation in the case of the stator copper temperature T s and cold air temperature T c a . The experimental data and the simulation are less well matched for the lower temperature region before t ≈ 300 s in the case of stator iron temperature T Fe and hot air temperature T h a . In the case of T Fe , the mismatch between the experiment and the simulation prior to t ≈ 300 s as shown in Fig. 9 is caused by the influence of the heat transfer parameter. The posterior distributions of the heat transfer parameter UA Fe2a and UA s2Fe related to stationary parts are more wider than the heat transfer parameter UA r2δ related to rotating parts. UA r2δ is the heat transfer related to the rotating copper, UA s2Fe is the heat transfer related to stationary copper, and UA Fe2a is the heat transfer related to stationary iron. From this result, we can infer the dispersion characteristics of the posterior distribution of the stationary copper, stationary iron, and rotating copper in the case of the air-cooled synchronous generator. The posterior distribution of the heat transfer parameters related to the stationary parts is more wider than the posterior distribution of the heat transfer parameters related to the rotating parts. The stationary iron receives heat from the heated stator copper at temperature T s because of terminal current I t flowing through the stator copper. Thus, the posterior distribution of heat transfer parameter related to stationary iron is more wider. This is because of the result from Section IV-D that the posterior distribution of the heat transfer parameter related to iron receiving heat from the stationary copper UA s2Fe and the heat transfer from heated stationary iron to air UA s2Fe are more wider. Thus, the simulation results of the stator iron T Fe is not matched with the experimental data for lower temperature region before t ≈ 300 s for the air-cooled synchronous generator. In the real operation of the air-cooled hydro generator at Åbjøra, Norway [15], for a period with the average time constant of 53 min the machine was in the state of the cold-run as described in Section II-B. The stationary iron takes time to get heated from the stationary copper heated from the terminal current I t .
Similarly, in the case of hot air temperature T h a , the mismatch between the experiment and the simulation before t ≈ 300 s as shown in Fig. 9 is because of the influence of the heat transfer parameter related to the stationary iron; UA Fe2a and UA s2Fe . From Fig. 1 the iron heat sourcė Q Fe2a = UA Fe2a (T Fe − T h a ) [16] which indicates that the nonhomogeneous heating of iron during the cold-run of the hydro generator cause heat transfer parameter related to stationary iron UA Fe2a to attain different modal values as shown in Fig. 5 (c). In addition, this non-homogeneous heating of the iron also affect the iron heat source parameterQ σ Fe as shown in Fig. 5 (e). This means that during the cold-run of the air-cooled synchronous generator the air inside the hydro generator at temperature T h a is heated intermittently or non-homogeneously as indicated by the experimental data during the cold-run prior to t ≈ 300 s. Similar, intermittency can be seen in the case of experimental data for cold air temperature T c a and stationary iron temperature T Fe . As the air temperatures T c a , T δ a and T h a inside the machine are interrelated through governing Eqs. (3), (4), and (5), it can be predicted that the experimental data during the cold-run of the hydro generator for air-gap temperature T δ a should also be intermittent. The governing equations for the metal and air temperatures formulating a dynamic model for the air-cooled machines with stationary iron parts need intermittent correction for its heat transfer parameter related to iron during the cold-run of the machine.

V. CONCLUSION
For the air-cooled synchronous generator as described in Section II-A, results from the analysis of the estimated parameters as described in Section IV-D show that the posterior distribution of temperatures of the stationary parts inside the air-cooled synchronous generator is more wider than the rotating parts. In the case of the heat transfer parameters, the posterior distribution of the heat transfer parameters related to metals is more wider than the posterior distribution of the airrelated parameters. Furthermore results also indicate that the posterior distribution of the heat sources parameters related to iron is also more dispersed than other heat sources parameters like heat source due to friction, etc.
From Section IV-E, results indicate that all the parameters estimated are practically identifiable. From Section IV-F, results indicate that the mismatch between the experimental data and the simulation results for the iron temperature and hot air temperature during the cold-run of the hydro generator is due to the higher dispersion characteristics in the posterior distribution of the stationary parts of the generator. The stationary iron takes time to get heated from the stationary copper. The stationary copper gets heated from the terminal current flowing through the stator copper. Thus, the posterior distribution of parameters for the air-cooled synchronous generator affects the mismatch between experimental data and the simulation results. Results also indicate that the heat transfer parameter related to iron attains an intermittent value during the cold-run of the air-cooled synchronous generator.

VI. FUTURE WORK
The governing equations for the air-cooled hydro generator, represented by Eqs. (1)(2)(3)(4)(5)(6), are considered with constant metal resistances and constant specific heat capacities. Future work includes parameter estimation and identifiability in the case of temperature dependent resistances and specific heat capacities. From Figs. 5 and 6, we see extremely narrower parameter distribution. It would be interesting to see the relationship between the measurement data available for the descriptor and the width of parameter distribution.
MADHUSUDHAN PANDEY received the M.Sc. degree in electrical power engineering from the University of South-Eastern Norway (USN), Norway, in 2019. Since 2019, he has been a Ph.D. Researcher with TMCC, USN, Porsgrunn, Norway. His research interests include modeling, control, and optimization of dynamic systems and integration of dispatchable renewable energy sources with variable renewable energy sources.
BERNT LIE (Member, IEEE) received the Ph.D. degree in engineering cybernetics from the Norwegian University of Science and Technology (NTNU), Norway, in 1990. From 1987 to 1991, he was an Assistant and an Associate Professor at NTNU. In 1992, he joined at the University of South-Eastern Norway, where he is currently a Professor in informatics and the Leader of the Telemark Modeling and Control Center. His research interests include control relevant modeling and advanced control, with applications mainly within the process and energy industries. VOLUME 10, 2022