HYDRO GENERATOR THERMAL SYSTEM IDENTIFICATION AND PARAMETER ESTIMATION

Original scientific paper In this paper identification and parameter estimation of a hydro generator thermal system using collected input and output data is presented. Third order state space model is constituted on the basis of a set of independent output variable (temperature) measurements derived from the system response to three thermal sources. Data were collected during a heat run and site acceptance test of a refurbished machine. Least squares (LS or RLS) methodology is used in parameter estimation procedure. Minimisation is performed on weighted measurement and model output distance squares. During the procedure system and input matrices were obtained. Thereafter, spectral properties of system matrix are analysed using calculated eigenvalues and eigenvectors. Eigenvalues assignment to the physical system objects is performed via L2 norm calculation and ranked to output variable sensitivities with respect to their time constants (negative reciprocal eigenvalues). Finally, simplified analogous lumped parameter scheme of the machine thermal system has been outlined.


Introduction
The main task of generator temperature monitoring system is to decrease temperature variations on critical locations regardless the actual electrical and magnetic load.It is important to reduce thermal dilatation of electric conductor insulation, especially out of core winding ends and thus to maintain the machine designed lifetime.Unfavourable temperature conditions, i.e. high temperature and its variations may lead to chemical and physical degradation of machine insulation and lifetime shortening.Design of such a machine temperature monitoring system requires very detailed knowledge of dynamic parameters and properties of the generator thermal system.Fundamentals of thermal static and dynamic processes can be found in [1,2].
Hydro generating unit maintenance scheduling is explained in [3].Detailed image consideration and investigations of temperature fields founded on physical laws with very developed geometric structure of a machine thermal model are given in [4÷7].Methodologies found in these sources relate to machine construction, thermal optimisation and temperature field distribution foremost by means of numerical approaches.Numerical model of a heat transfer in a finned housed electrical motor is validated experimentally using thermography in [8].Example of elaboration for thermal steady state of an asynchronous machine is given in [9].In this reference second order system model of temperature dynamics in asynchronous machine is exploited.Transient thermal phenomena are coupled with two main machine construction corpora: copper electric conductors and iron core as magnetic conductor.Generally, the dominant thermal dynamics in all electric machines relates to these corpora.Similarly, detailed model of total enclosed fan cooled induction motor thermal behaviour is given in [10].The cooler structure of a high voltage asynchronous motor is optimized using the theory of Computational Fluid Dynamics in [11].In special cases, when machines are supplied from variable frequency sources, frequency influence on thermal processes must be carefully taken into account [12].
Since analytical modelling of electrical machine thermal system is a very complex issue, in this paper the new hydro generator thermal model based on parameter identification and estimation is presented.Referent state space identification and parameter estimation procedures of linear signals and systems (ARA, ARMAX) based on their detailed observations, especially of discrete systems can be found in [13] and [14].Based on such methodology the authors in [15] give lumped parameter estimation for fourth order thermal system of a semiconductor in order to predict the temperature dynamics in high current load conditions.They applied recursive parameter estimation procedure.
A least squares methodology for induction motor parameter estimation is used in [16].Similar methodology is also used in this paper.After explaining parameter estimation procedure, hydro generator thermal system model is built on the basis of independent temperature measurements derived from the system response to three thermal sources.Measurements are performed during a Technical Gazette 24, Suppl.2(2017), 257-264 heat run and site acceptance test of a refurbished machine.Model quality is examined by comparing estimated and measured temperatures.Besides avoiding complex modelling of system behaviour, the main advantage of proposed method is possibility to assign gained eigenvalues to appropriate physical system objects via L 2 norm ranking of system time response sensitivities with respect to parameters, i.e. time constants.

Thermal model and parameter estimation approach
Parameter estimation procedure for linear discretized system is based on general discretized state space system model given in [13,14,15] and [17]: ).
Θ Θ (6) Relatively complex estimation procedures, even for linear system generally exposed to disturbances which may have very different frequency spectra that must be taken into consideration, especially if input and output measurement sets are extremly large, are explained in [13,14,15] and [17].These frequency spectra should not have consideralbe impact on estimated parameter accuracy and bias.
Simplified expression of parameter matrix Θ estimation, according to Eq. ( 5), is given in the following way: F is regression matrix that is formed from output vector y d and input vector u d , taking into account definition of linear discrete form ( 1) and ( 2).After calculation process ( 7) is finished, matrices Φ and Γ are found by extraction from Θ as given in Eq. ( 5).The model response y ˆ is calculated according to ( 1) and (2) using mesured input vector u d : ) Their continous equivalents (A c and B c ) will be found for known sample interval T s .Doing so, it will be possible to calculate eigenvalues of continous system matrix A c .The other approach to calculate eigenvalues of c A consists of calculating eigenvalues of Φ first and recalculating them directly for continuous equivalent using the same sample interval T s .If gained eigenvalues are real and negative numbers, then their negative multiplicative inverse values are corresponding to some time constants.After the Eq. ( 6) is successfully minimised and matrices (A c and B c ) calculated, the differential equations that satisfactorily interpret real system dynamics can be written: . d Predominant belonging of a particular eigenvalue to constitutive part of a system can be determined via L 2 norm ranking of calculated system response sensitivities with respect to above mentioned time constants in the next manner.Subscript of measured or model variable (y d or y) precisely indicates the system constitutive object (e.g.copper temperature means copper as object).According to Eq. ( 9) it is possible to obtain time output vector sensitivity with respect to T i .If the system has n eigenvalues, i.e. time constants, then the sensitivity system response contains n 2 scalar time functions or: .
For example, if the response y k (t) has the largest L 2 sensitivity norm (s max ) with respect to T i : ( ) (11) then time constant T i predominantly belongs to the y k (t), i.e. to the k th object described with this variable y k (t).
In other words, parameter T i predominantly belongs to the variable having the largest time response change when the parameter T i changes for some finite difference ΔT i : max ( , ) Namely, model description according to Eq. ( 9) can be achieved by considering physical phenomena of thermal processes represented by analogous scheme with lumped parameters.Basically, these parameters are thermal resistances (R θ , K/W) and thermal capacitances (C θ , Ws/K).More details about this representation will be given later in paper.

Measured values type and quantity
Machine heat run test starts at synchronous speed with its excitation, synchronisation to the utility transmission network and finally with achieving rated load.Rated load must be maintained until the temperature steady state is practically attained.Duration of the loading state should not be shorter than five times the largest time constant.In our case the largest time constant corresponds to the iron core thermal time constant (to be determined later).Specifically, thermal system observation lasting was 20868 seconds (about 5,8 hours) with 10 ms sample interval, but for analysis the recordings are decimated 100 times, so the sample interval was enlarged to value T s = 1 s.Fastest phenomena would allow larger sample interval, even greater than 5 s.Hydro generator cooling system is shown in Fig. 1.Cold water supplied from a cold water tank is brought to generator heat exchanger.Cold air from the heat exchanger is used to cool the turbine and generator bearings, as well as the generator windings.Measured quantities recorded in digital form and saved for further analysis are divided into output and input variable groups.
Group of output variables, that form thermal system state variable set representing response to the certain input variables set applied to the system, consists of:  Certain objection to this set of variables can be addressed to rotor winding temperature measurement absence.Rotor winding temperature would also constitute one independent state variable.This temperature can be obtained via excitation control system, as it can be estimated on the basis of rotor winding resistance change with temperature.This quantity will be taken into account in the future measurements and estimations.
Second group of measured quantities contains the set of variables used for input variables calculations: Following quantities are also measured: • Cooling water flow, Q, l/min,

•
Cooling system inlet cooling water temperature ("cold water"), θ ccw ,°C.While all temperatures are measured via PT100 probes, voltage, currents and water flow are collected from corresponding transducers.Input variables are formulated as heat sources applied to the machine thermal system.First variable of input vector u d is calculated using the armature current and represents the stator winding copper losses p Cu_a in the form: It is important to note that the temperature change of armature resistance was taken into account.Second input variable is obtained in a similar manner and represents field winding copper losses p Cu_f : Although rotor winding temperature θ fCu is not directly observed, its final temperature is known to be 81 °C.Knowing this, rotor winding temperature rise can be obtained from armature winding temperature rise multiplied by a factor that ensures the final temperature value equal to 81 °C.Technical Gazette 24, Suppl.2(2017), 257-264 Third input variable represents iron core losses using generator voltage instead of magnetic induction or flux linkage.From rated iron losses (Tab. 1) equivalent resistance is obtained and voltage quadratic law was assumed so that p Fe can be calculated: W. , Additional copper and iron core losses (Tab. 1) for rated operation point were divided by two and consistently added to Eqs. ( 14), ( 15) and ( 16), respectively.Input variables ( 14), ( 15) and ( 16) can be expressed in pu values dividing them by apparent rated power (S n , VA, Tab. 1).
Power leaving the thermal process at the end of observation dW cw /dt, expressed in Watts, is proportional to output cooling water final temperature rise (Δθ cw = 3,1 °C) multiplied by its flow (Q) and is calculated in the following way: Measured water flow Q (l/s) was rather constant and its average value was about 48,2 l/s.Quantity, c = 4,183×10 3 J/(lgK) is specific water heat at 25 °C, and ρ = 1 kg/l is water mass density.The heat dissipated through the cooling system (17) can also be given in pu values dividing it by generator rated apparent power (S n , VA, Tab. Figure 2Armature and field winding copper losses, iron core losses and cooling system as excitation applied to hydro-generator thermal system Eigenvalues of equivalent continuous system matrix according to [17]  Time behaviour of three input variables as thermal system excitation is shown in Fig. 2. Time responses of real system and system model are presented in Figs. 3, 4  and 5, respectively: temperature rise of armature (copper) winding (Δθ Cu , pu), iron core temperature rise (Δθ Fe , pu) and output cooling water temperature rise (Δθ cw , pu).As it is shown, 1 pu corresponds to 80 °C.
Values relate to modelling quality of copper, iron core and cooling water temperature rises respectively.These values and time dynamics in Figs. 3, 4 and 5 indicate an excellent agreement between measured and simulated system response.
Calculated eigenvalues of continuous system (23) belong to the entire system, not to particular constitutive object.However, based on consideration in chapter 2, it is possible to determine predominant belonging of particular eigenvalue to certain system object or the state variable.Instead of eigenvalues it is more convenient to use time constants (in seconds) as negative inverse of eigenvalues: Each response is at most affected by particular time constant change.This influence can be quantified in sense of already introduced Eq. ( 11).Fig. 6 shows that the parameter T 1 = 23,901 s has the highest influence on cooling water temperature rise (variable 3 according to the state variable definition in Eq. ( 13)): Consequently, time constant T 2 = 737,96 s mostly belongs to the armature winding, i.e. copper.
Finally, time constant T 3 = 3344,5 s predominatly belongs to the iron core temperature rise (Fig. 8, Fe, iron core, variable 2 according to the state variable definition in (13)).Using the system matrix A c , Eq. ( 22) and system input matrix B c , Eq. ( 23), under assumption that system matrix A c is not singular and steady state condition of Eq. ( 9) is fulfilled, final values of state variables y(∞) can be calculated: . ) ( 0 d  This calculation shows that thermal system is sufficiently long observed in the sense of achieving practical steady state temperature conditions.Similar results were obtained using instrumental variable approach (IV method).
We can return to the structure of matrices A c and B c , given in Eqs. ( 22) and ( 23).All elements are different from zero.Obtained system is of third order with three input variables.This structure approximately corresponds to analogous lumped parameter scheme of hydro generator thermal system represented in Fig. 9.
Also, distributed thermal sources and sink are shown in this scheme via k ij coefficients, such that i = 1, 2, 3, j = 1, 2, 3. Relation between these coefficients and B c matrix coefficients could be established by writing differential equations directly from Fig. 9. On the basis of the least square (LS, or RLS) methodology, hydro generator thermal state space system model has been identified and parameters were estimated.Estimated parameters satisfactorily interpret the most substantial spectral properties of the real system.Records of digitally measured and recorded collection of ten relevant variables during 5,8 hours observation of machine heat run and site acceptance test were available.Final set of measured data is reduced to six independent process variables.Three of them describe temperature dynamics of armature winding, iron core and output cooling water, whilst remaining three variables represent system thermal excitation.
Predominant belonging of particular spectral parameters (eigenvalues, i.e. time constants) to certain constitutive part of a system was based on L 2 norm ranking of system time response sensitivities with respect to parameters (time constants).In this way the dominant system time constants were determined.Obtained system model enables prediction of temperature behaviour for machine constitutive parts for different electric and magnetic loads.Often such measurements cannot be obtained at site.Example of such measurements is very accurate verification of machine rated apparent power under referent thermal conditions related to maximal temperature constraints deduced from temperature standards for thermal and insulation class in electric machines manufacture.Thus the model obtained can be used for machine temperature monitoring, system design and corresponding temperature controller tuning.
Further improvements of model by means of introducing field winding temperature as an addition to estimation model are possible.Further research will be focused on accurate generator inner space temperature observation.

Figure 1
Figure 1 Hydro Generator Cooling System: CW -cold water, HW -hot water, CA -cold air, HA -hot air

Figure 3
Figure 3 Measured and calculated responses of armature winding (copper) temperature rise

Figure 4 Figure 5
Figure 4 Measured and calculated responses of iron core temperature rise

Figs. 6 ,
Figs. 6, 7 and 8 subsequently show the system response sensitivities for small changes in time constants T i , for constant T 1 = 23,901 s mostly belongs to the cooling system, i.e. the cooler.Further, parameter T 2 = 737,96 s has the highest influence on armature winding (Cu, copper) temperature rise (Fig.7, variable 1 according to the state variable definition in (13)): the system parameters are estimated successfully, all final values (t = ∞) of temperature rises can be calculated.In this case, it is sufficient to use input vector in the last time observation point (u d (k) = u d (20869)):

( 32 )Figure 6 Figure 7 Figure 8
Figure 6 System time sensitivity responses with respect to time constant 1 23,901 s T =

Figure 9
Figure 9 Analogous scheme of hydro-generator thermal system with lumped parameters: thermal resistances ( θ R ,°K/W) and thermal capacitances ( θ C ,Ws/°K) corresponding to Ac and Bc structure

Table 1
Hydro generator manufacturers data