Identification of Parameters for a Simplified Model of Large Rotating Electric Machines Using Time Series and Updating Techniques
S. Silva1, S. B. Shiki1, F. Franzoni1, G. C. Brito Jr.1,
2 1 Western Paraná State University (UNIOESTE), Centro de Engenharias e Ciências Exatas (CECE) Grupo de Pesquisa em Din?mica de Estruturas e Máquinas (GPDEM) Av. Tancredo Neves, n.? 6731, Parque Tecnológico Itaipu (PTI), Foz do Igua?u, Paraná, Brazil email: sam.silva13@gmail.com, sbshiki@gmail.com, felipe.franzoni2@gmail.com Itaipu Binacional Av. Tancredo Neves, n.? 6731, Foz do Igua?u, Paraná, Brazil email: gcbrito@i taipu.gov.br
2
Abstract
This paper describes the application of simplified analytical mathematical models of the dynamic behavior of large rotating electric machines from Itaipu Power Plant. The parameters in this model are updated based on real data obtained in the commissioning of generating units of Itaipu. Auto-regressive moving average models (ARMA) are identified with the output time series measured at the bearings stations. The modal parameters are estimated from these ARMA models and used as reference for updating the stiffness parameters in the simplified model of the hydro generating units. Genetic algorithms are implemented to find the optimal solution of the residuals between the differences from analytical and experimental modal parameters. Sensitivity analysis is performed to define the best parameters of stiffness to update in the analytical model. The results show that a representative and simplified model can be used for primary dynamic analysis of operation and future health condition monitoring.
1
Introduction
Itaipu Power Plant has 20 generating units where each one provides about 700MW. The dynamic analysis in the rotating electric machines in Itaipu roles a strategic function, as this power plant generates 90% of the electric energy consumed in Paraguay and more than 20% of the Brazilian electric grid demand [1]. So, the study of algorithms for damage detection, diagnosis and prediction in the Itaipu Power Plant is an area of interest, mainly due to the demand for better performance, reliability and availability requirements to electrical systems of both countries. The use of mathematical models to describe the dynamic behavior of the hidrogenerator is essential, since it would allow more control of the structural integrity of these machines. However, normally the mathematical models are not be able to set correctly and be compatible with experimental behavior [2]. In this sense, model updating techniques have been used to minimize the differences between analytical models and the experimental results. Several studies use different model updating technique, for example Pereira et al [3] proposes a methodology for structural damage detection based on a comparison between a reference finite element model and experimental frequency response function (FRF). If the system presented some structural changes, the modal parameters also changed and, consequently the FRFs modified. Kim and Park [4] discuss the need to approach the model updating process as a multiobjective optimization problem and proposes a genetic algorithm based on the ParetoOptimal concept on this type of problem. Zivanovic et al [5] describe the complete process of the finite
element modeling and updating of the Podgorica bridge. The same uses a sensitivity analysis of the dynamic responses to decide the best parameters to be used in the model updating process. In the present paper an update procedure in a simplified mathematical model with two degrees of freedom to describe the dynamic behavior of the generating unit 09A. An autoregressive moving average model (ARMA) is used to identify the natural frequencies of the machine through vibration data measured in the station bearings. The most sensitive model parameters are updated to minimize the differences between the natural frequencies of the analytical model and experimental natural frequencies. This inverse problem is solved with the use of sensitivity analysis and genetic algorithms. The results found are discussed in order to show the applicability of this procedure
2
Simplified model of generating unit
Brito Jr. et al [6] modeled only the 09A generating unit to analyze the dynamic behavior of this rotating machine. This model has two degrees of freedom, the axial displacement (x(t)) and the angular displacement (θ(t)) as showed in the fig. (1).
Figure 1: Detail of the degrees of freedom of the machine and his physical model. The mathematical model, disregarding the effect of damping can be written as:
? m11 ?m ? 21
&(t )? ? k11 m12 ? ? & x ? && ? + ? m22 ? ? ?θ (t )? ?k 21
k12 ? ? x(t )? ? F (t ) ? ? ?=? ? k 22 ? ? ?θ (t )? ?? Z FG F (t )?
m 22 = J = 1,36 × 10 8 kg ? m 2
(1)
The inertial parameters in eq. (1) are:
m11 = mTotal = 2780ton m12 = m 21 = 0
where mTotal is the total mass of the machine, composed by the mass of turbine rotor including the water and mass of generating rotor, and J is the moment of inertia relative to the center of mass. The stiffness parameters in eq. 1 are written as: k11 = k 0 + k1 + k 2 + k 3 k12 = ?k 0 Z 0 g ? k1 Z 1g ? k 2 Z 2 g ? k 3 Z 3 g k 21 = ? k 0 Z 0 g ? k1 Z 1g ? k 2 Z 2 g ? k 3 Z 3 g
2 2 2 2 2 k 22 = k 0 Z 0 g + k1 Z 1 g + k 2 Z 2 g ? k 3 Z 3 g + k 5 R
9 9
(2) (3) (4) (5)
where k 0 = 0.83 × 10 N/m is the radial stiffness of generator upper guide bearing; k1 = ?0.6 × 10 N/m is the generator radial magnetic stiffness; k 2 = 1.25 × 10 9 N/m is the radial stiffness of generator lower guide bearing; k 3 = 1.43 × 10 9 N/m is the radial stiffness of turbine guide bearing; k 5 = 32.0 × 10 9 N/m is the thrust bearing axial stiffness, z 0 g = 7.8 m, z1 g = 4.8 m, z 2 g = 2.8 m, z 3 g = ?5.7 m are the distances relatives to
center of mass, and R = 2.1125 m is the radius in the thrust bearing pads pivoting point . These values were obtained by tests and project values so that are subject to uncertainties. Although simplified, this model showed to be useful for a first analysis and the monitoring of the conditions of these rotating machines. Silva and Brito Jr [7] showed a damage detection method based on the ARMA model and statistical process control (SPC), using the simplified mathematical model to simulate the radial vibrations.
3
ARMA model
Each set of samples from the radial vibration of the bearing gi(k), i = 1,2, K , m and k = 1,2, K , n , represents the time series for m positions sizes in n discrete instants of time. Each time series gi(k) is normalized to remove trends [8]: g (k ) ? ? i (k ) (6) xi (k ) = i σ (g i )
where xi(k) is the ith signal with zero mean and unitary standard deviation, gi(k) is the sampled signal and ?(gi) and σ(gi) are the mean and standard deviation operators, respectively. To improve the frequency resolution in a low band, a decimation procedure is normally required [9]. An auto-regressive moving average model (ARMA) can be described as [10]: A q ?1 xi (k ) = C q ?1 ei (k )
th
( )
( )
(7)
where e(k) is the i error between the measured signal and the output from ARMA prediction model and q ?1 is the delay operator in the polynomials. A(q-1) and C(q-1) are written by the following equations:
( ) C (q ) = 1 + c q
?1 1
A q ?1 = 1 + a1 q ?1 + L + a na q ? na
?1
(8)
+ L + c nc q
? nc
where na and nc are the delay orders in the polynomials A(q-1) and C(q-1). The choice of this orders are essential for the fit and prediction of the ARMA model. In this paper the Akaike information criteria (AIC) was used to find the best orders for the model [11]. The prediction error method (PEM) can be used for estimating the coefficients A(q-1) and C(q-1). This method is already implement in the system identification toolbox for Matlab? and Scilab?. To validate the model two methods can be used, the first consist in a direct comparison between the model results and experimental data. The second is a residual analysis with the cross-correlation function (CCF) between the residues of the model and the output signal. In both validation procedures, several data in different days are used. The roots of the polynomial A(q-1) are directly related with the natural frequencies and damping factors of a linear system, However, these poles are in the z-domain (discrete model). The discrete poles in the z-domain must be converted to s-domain (Laplace) through a mapping procedure, as for instance the bilinear transformation [9] to be correlated with the modals parameters. The relationship between the continuous poles and the modals parameters (natural frequency and damping factor) can be described by: (9) s = ?ζ ω ± ω 1 ? ζ 2 j
n n n n n
where sn are the continuous poles of the transfer function, ζn the damping factors, ωn the natural frequencies and j is the complex number.
4
Model updating
A problem often encountered in engineering is called inverse problem, where the model of a system must be constructed or refined given some input and/or output data from the system [2]. Many engineering activities can be considered inverse problems: model updating, experimental modal analysis,
parameter estimation, damage identification, input forces identification, etc. Model updating seeks to correct differences between the analytical model and the dynamic responses measured in vibration tests done in a structure. The general procedure of this method is showed in the flowchart in fig. (2).
Analytical model Yes Correlated? Updated model
Experimental data No
Model updating
Figure 2: Flowchart of the model updating process.
The differences encountered in the analytical model can be originated from various sources [12][13]: Simplifications of the theory adopted to derive the model; Difficulty in modeling dissipative effects (damping); Discretization errors; Numerical errors occurred during the resolution of the equations of motion; Uncertainties of the values of some parameters of the system used as inputs to generate the model, resulting from the difficulty in estimating accurately the values of physical and geometrical properties. The methods of model updating can be classified in two classes: the direct methods and the iterative or parametric methods [5]. The first class is based on updating of stiffness and mass matrices directly, in a way that often has no physical meaning. The second class of iterative methods, on the other hand, concentrates on the direct updating of physical parameters which indirectly update the stiffness and mass matrices in which these parameters are inserted. Although computationally more expensive, the iterative methods are more used because its characteristics lead to corrections in the model that are easier to be interpreted.
? ? ? ? ?
5
Sensitivity analysis
The mass, damping and stiffness matrices have all the information of the physical and geometrical characteristics of the modeled structure [12][14]. However, an important information for the model updating process is the degree of the variation of the dynamic behavior as a result of changes in the physical and geometrical properties. Hence, the sensitivity analysis aims to establish the relationship between the dynamic responses and the variation of certain structural parameters. In the mathematical point of view, the process of this analysis corresponds to compute the partial derivatives of the model responses with respect to each parameter of interest. Thus, the sensitivity of the dynamic response r with respect to a certain parameter pi, to a value pi=pi0 is defined by:
?r ?p i
= lim ?pi →0
pi0
r M ( p i0 + ?p i ), C ( p i0 + ?p i ), K ( p i0 + ?p i ) ? r M ( p i0 ), C ( p i0 ), K ( p i0 ) ?p i
(
) (
)
(10)
where M is the mass matrix, C is the damping matrix, K is the stiffness matrix and ?pi is a small increment in the value of pi0. This value can be computed in two ways: by computing the analytical partial derivatives that can be sometimes computationally expensive, or calculating the derivatives by using the finite differences approach. Zivanovic et al [5] and Bakir et al [15] used sensitivity analysis to help in the choice of the parameters to be updated in the model updating process. Parameters with null sensitivity or next to zero
are discarded because for these it would be necessary great modifications to obtain minor changes in the structure dynamic responses. So, it is better to choose the values with high sensitivity (positive or negative) since its effect are more intense in the model responses.
6
Inverse problem solution
Optimization problems can be defined as the determination of variables so that an objective function reaches an extreme value (maximum or minimum) subjected to some constraints of the problem. Algorithms to solve this kind of problems are often classified in two groups: the classic methods based on the calculation of gradient values (derivatives) that provide the search direction of the algorithm [16], and the heuristic methods that changes the optimization parameters based on random decisions [17]. Despite the popularity of the classical methods, often is not possible to ensure that the final solution found by these strategies is actually the global optimum. That characteristic is dependent on the level of complexity of the optimization problem [18]. In these cases is usual to use heuristic algorithms as: genetic algorithms, simulated annealing, particle swarm, etc. Genetic algorithms are methods of search and optimization that simulates the natural process of evolution, by means of the natural selection of the species described by Charles Darwin [19], [20]. These algorithms are robust methods and applicable in the resolution of various problems. In the general procedure of a basic genetic algorithm, in the first step, an initial population of chromosomes or individuals is created which represents possible solutions to the problem codified usually in binary code. In the next step, each chromosome is evaluated by a quality measure called fitness that is related to the objective function of the problem. After this, crossover and mutation operators are applied generating a new population of chromosomes. This process is iteratively repeated until a pre-established stop criterion is reached or until a maximum number of generations. The flowchart of the basic genetic algorithm is showed in fig. (3).
No
Initial population Evaluate fitness
Reached the stop criteria?
Crossover
Mutation
Yes
Best chromosome is the final solution
Figure 3: Flowchart of a basic genetic algorithm.
6.1
Multiobjective optimization
In most practical situations the optimization problems have several objective functions that may be minimized or maximized simultaneously. This characteristic brings some difficulties since this evaluation functions can be conflicting. The study of multiobjective optimization is important because the model updating problem often fulfills these criteria. Many objective functions can exist based in the natural frequencies, mode shapes and FRFs [4]. There are many techniques to deal with the multiobjective problem [12], but in this work it will be described two techniques often founded in the literature. The first methodology consists of combining the objective functions with a weighted sum [5]:
1 2 2 ?K ? ? wk × f k ( x ) ? f k' ( x ) ? ? ? f (x ) = ? ? ? ? ' ? f pv ( x ) ? f k ( x ) ? ? ? k =1 ? ? ? ?
∑
(
)
(11)
where f(x) is the resulting objective function, wk is a weight factor, fk is one of the objective functions to combine, f’k is the desired value of the objective function and fpv is the worst value accepted to the objective function. The second methodology uses the concept of Pareto-Optimal Set [4], [19], [20]. This concept can be explained by a simple example. Consider a minimization problem with two different objective functions F1 and F2, and five points to be evaluated: A, B, C, D and E with the functions in the two coordinates as showed in fig. (4).
Figure 4: Illustration of the concept of Pareto-Optimal set.
Scanning the graph further reveals that the best points are next to the origin (0,0). It may be noted that in the group of the points A, B and C none of the three points is best along both dimensions, there is gain along one dimension and loss along others. These points are called nondominated because there are no points better than these on every criterion. The remaining points have worse values than the first group and are called dominated. This method consists in obtain an optimum group of nondominated values called Pareto-Optimal Set instead a unique solution.
7
7.1
Results
Estimation of modal parameter by using experimental data
The signals of vibration of the generating unit 09A (Itaipu) are measured in three different positions: at the generator upper guide bearing, at the generator lower guide bearing and at turbine guide bearing. Each point has two sensors placed in orthogonal localization, called X(Brazil) and Y(upstream), with six signals of vibrations fig. (5). The data are relatives to the radial vibration of the shaft measured in ?m. The sampling rate was 1,2 kHz with 29977 samples for each sensor installed. The tests were conducted under permanent conditions. The signals are processed by an analog low pass filter with cutoff frequency of 1 kHz.
Figure 5: Unit 09A with the measured points.
Before identifying ARMA models, all signals are normalized to remove trendy. A downsampling procedure is also realized in order to reduce the maximum frequency (N = 2998 samples). After these procedures, the data are divided in two parts, the first is used to construct ARMA models and the second one to validate. With the data properly sampled ARMA models are constructed from each set of samples. The polynomials were estimated by the prediction error method [10] and the orders regression were chosen with Akaike information criteria [11], which was used interactively for each set of samples testing a large number of orders. These results were analyzed to find the order that well represents all the signals. The polynomials A(q-1) obtained are characteristic equation of transfer functions. In this paper was considered a mean value of each identified frequency to determinate their origins, which can be from the generating unit or external factors. Some physical characteristics from generating unit (60 Hz in Brazil) are used to eliminate known harmonics with no relatives to model behavior, see tab. (1): ? ? ? ? the nominal rotation (92,3 RPM); the number of blades of the turbine (13 blades); the number of guideline blades of the distributor (24 blades); the number of poles of the generator (78 poles);
Origin of the frequency Rotation frequency and its harmonics Electromagnetic frequencies Frequency of blade passage Frequency of guideline blade passage Frequency of the generator poles passage f [Hz] 1,54 60, 120, 240 20 36,92 120
Table 1: Generating unit frequencies from physical characteristics. Figures (6), (7), (8) and (9) show the described procedures executed for one of the signals (upper guide bearing in Y direction). Figure (6) show the sampled signal after normalized and downsamplig procedure. Figure (7) show the interactively Akaike procedure which provided an minimum index for 14 and 6, na and nc, respectively. The direct comparison shows a good fit (95,51%) between the prediction output and the measured output, see fig. (8). The residual analyses illustrates that the ARMA model can represent well the dynamics contained in the data, because the correlation is into a range of tolerance. So, this confirms that the model residues have random features, see fig. (9).
Normalized sampled signal
4 3
?5
AIC
Radial displacement ? x(t) [?m]
2
?10
1 0 ?1 ?2 ?3 0
AIC index
?15 ?20 ?25 ?30 0 5 10 15 10 20 15 5 0
5
10 Time [s]
15
20
25 A(q) order
C(q) order
Figure 6: Normalized sampled signal.
Figure 7: Akaike information criteria.
Direct comparison
1.5 datavalid; measured model; fit: 95.51% 1.2 1 0.8 0.5 0.6 0.4 0.2 0 ?1 2 4 6 Time [s] 8 10 12 ?0.2 0 5
Cross correlation function
Radial displacement ? x(t) [?m]
1
0
?0.5
CCF
10
15
20 Lag [N]
25
30
35
40
Figure 8: Direct comparing between the sampled signal and the model signal.
Figure 9: Cross-correlation between the sampled signal and the model residue.
Table (2) shows the 1st natural frequency to all signals identified with experimental data:
Signal Upper bearing of the generator in X Upper bearing of the generator in Y Lower guide bearing of the generator in X Lower guide bearing of the generator in Y Guide bearing of the turbine X Guide bearing of the turbine Y Mean f1 [Hz] 4,62 4,55 4,54 4,62 4,56 4,55 4,57
Table 2: 1st natural frequency identified for each signal. The mean value is used as exact value and will be used to update the 1st natural frequency of the mathematical model.
7.2
Estimation of modal parameter using simulated signals
To simulate the signals of vibration with the mathematical model, the state space representation was chosen and the motion equations were solved by the state transition matrix method [21]. Based on [11] the sampling rate was considered as 30 times greater than the maximum frequency of interest, around 199 Hz, the number of samples was chosen as 4096. The parameters in the mathematical model in this simulation are the nominal design values, subjected to uncertainties. The input F (t ) is chosen based on the experimental displacement measurements in the machines. The order of this value is 106 ~ 107 N and the excited frequencies were only considered the fundamental frequency and ? of the frequency (excitation caused by partial load vortexes):
? ? 1 ?? (12) F (t ) = F ? ? sin (?t ) + sin? 4 ?t ? ? ? ? ?? ? where the frequency to the machine operating in 60 Hz (Brazilian electrical system) is 92,3 RPM. Four different data were simulated by this mathematical model in this work, where two represents the radial displacement and two represents the angular displacement. For all cases were added 10% RMS noise in the output signals and also all the signals were filtered by IIR band pass filters. The operational modal analysis was performed using discrete models (ARMA) constructed and validated based on simulated signals. The procedures to determine the ARMA model and correlate the poles of this model with the modal parameters were the same as the previous item of this paper. In this case a stabilization criterion to determine the frequencies from physical poles of the system was used [22]: it was considerer that the natural eigenfrequencies of a mode of the current and previous models have an absolute and maximum deviation of 1%.
Figure (10) shows the filtered simulated signal and the fig. (11) shows the implemented Akaike information criteria. Figures (12) and (13) show the validation procedures for one of the ARMA models with na 6 and nc 4 with direct comparison and cross-correlation function, respectively. The results were satisfactory, once the fit of the model is within a range considered acceptable and the cross-correlation function shows the random nature of the prediction error.
x 10 2
?4
Simulated signal
AIC
Radial displacement [?m]
?22.5
1
?23
AIC index
?23.5 ?24 ?24.5
0
?1
?25 1 2 2 3 4 5 3 6 4 1
?2 0 2 4 Time [s] 6 8 10
C(q) order
A(q) order
Figure 10: Sampled signal.
x 10 2.5
?4
Figure 11: Akaike information criteria.
Cross correlation function
1.2 1 0.8 0.6 0.4 0.2 0 ?0.2 0
Direct comparison
datVALID; measured ARMA; fit: 96.37%
Radial displacement ? x(t) [?m]
2 1.5 1 0.5 0 ?0.5 ?1 ?1.5 ?2 1 2 3 4 5 6 Time [s] 7 8 9 5 10 Lag [N] 15 20
Figure 12: Direct comparing between the sampled signal and the model signal.
For each signal (radial and angular displacements) just one eigenfrequency respected the stability criteria. Table (3) shows the value for these two frequencies for each model that was simulated.
Discrete model ARMA(na,nc) ARMA(4,2) ARMA(4,3) ARMA(4,4) ARMA(5,2) ARMA(5,3) ARMA(5,4) ARMA(6,2) ARMA(6,3) ARMA(6,4) Mean Maximum deviation [%] Radial displecement signal x(t) ω1 [Hz] 5.15 5.14 5.13 5.15 5.14 5.14 5.14 5.14 5.14 5.141 0.214 Angular displacement signal θ(t) ω2 [Hz] 6.64 6.63 6.64 6.64 6.63 6.65 6.63 6.63 6.63 6.636 0.211
Table 3: Eigenfrequencies identified from the signals of the mathematical model.
CCF
Figure 13: Cross-correlation between the sampled signal and the model residue.
Figure (14) shows the bode diagram for two ARMA models, one constructed with radial displacement signal and the other with angular displacement signal both with na = 6 and nc = 4.
Frequency response function (FRF)
From: v@y1 To: y1
20 0 ?20
System: Hcy I/O: v@y1 to y1 Frequency (Hz): 5.15 Magnitude (dB): ?39.2
Angular displacement Radial displacement
Magnitude (dB)
?40 ?60 ?80 ?100 ?120 Frequency (Hz)
System: Hc I/O: v@y1 to y1 Frequency (Hz): 6.66 Magnitude (dB): ?18.8
10
1
Figure 14: Bode diagram for radial and angular displacement ARMA models.
The mean value of each eigenfrequency is taken as the exact value of these parameters and they are used in optimization process.
7.3
Sensitivity analysis
In order to implement the sensitivity analysis it was employed the eq. (13) of normalized sensitivity proposed by [5] which calculates this value using the progressive finite difference approach: S norm
i, j
=
?Ri Pj ?Pj Ri
(13)
where ?Ri is the variation of the dynamic response, ?Pj is the variation of the parameter, Pj is the value of the parameter to calculate the sensitivity and Ri is the respective response. The equations 13 allow the comparison of the sensitivity values of different parameters. In this paper the sensitivities of the stiffness values (k0, k1, k2, k3 and k5) showed in eq. (2) to eq. (5) were compared, the dynamic responses used were the first and second natural frequencies. The variations of the responses were calculated with an increment of 1% in the value of each of the stiffness parameters. Figure 15 shows the sensitivities for each parameter and each dynamic response.
0.3 0.25
Normalized sensitivity
0.2 0.15 0.1 0.05 k1 0 ?0.05 ?0.1 ?0.15 1 2 3 First natural frequency Second natural frequency 4 5 k0 k2 k3 k5
Stiffness values
Figure 15: Bargraph showing the sensitivities of the natural frequencies to the stiffness values.
By observing the fig. (15), it could be concluded that the first natural frequency is more sensitive to variations in k2 and the second natural frequency is more sensitive to the parameter k5.
7.4
Updating of the simplified model of the generating unit
The procedure of model updating of the simplified model described in Section 2 was done by using a genetic algorithm. To update the stiffness parameters it was used the following equation:
k nadj = k nor × (1 + α n )
(14)
where knadj is the adjusted value of the n stiffness parameter for n ∈ {0,1,2,3,5} as described in Section 2, knor is the original value of the stiffness parameter n of the unchanged vibration model and αn represents the percentage which the original stiffness value increase or decrease. The genetic algorithm codifies the α values in binary vectors that represents real numbers with floating points, which are the chromosomes or individuals of the algorithm. The objective functions of this inverse problem are based on the relative error of the natural frequencies. In this work, two studies were made to show the capacity of the model updating.
7.4.1
First study: single objective optimization with experimental data
In this study, it was used the first natural frequency identified with the experimental vibration data in Section 7.1. This case has only one objective function because the signal used to identify the frequencies has only information of the first natural frequency. So, the model updating problem using this experimental data becomes a single objective optimization. For the first natural frequency the sensitivity analysis, showed that the first two best parameters to be updated are k2 and k3. These two parameters are updated to minimize the error between the first natural frequency obtained with the experimental signal and the frequency of the updated model analytically obtained with the resolution of a problem of eigenvalue and eigenvector. The objective function can be written as:
Min J ( p (α ) ) = ω1adj ( p (α ) ) ? ω1 exp
, for -50% < α < 50%
(15)
ω1 exp
by the solution of a problem of eigenvalue and where ω1adj is the natural frequency obtained eigenvector based on eq. (1) with the stiffness values modified, ω1exp is the frequency obtained with the processing of the experimental signal, p is the vector with the k2 and k3 stiffness values that will be updated, α is the vector with the α2 and α3 correction values which the range are set in -50% to 50% to ensure a small modification of the parameters and J is the objective function or the fitness of chromosomes. In this optimization problem the genetic algorithm codifies α2 and α3 values in binary vectors that compose the chromosomes. The genetic algorithm used to solve the inverse problem has the configurations showed in Table 4.
Type of selection Type of crossover Range of α Digits of precision Population size Number of participants of the tournament Number of iterations (generations) Crossover tax Mutation tax Tournament One point [-0.5 0.5] 3 100 5 100 0,8 0,05
Table 4: Genetic algorithm configurations In order to ensure that the fitness of best chromosome in each generation has been selected, an elitist procedure was implemented [20]. In the present study, the best chromosomes founded in 10 different simulations with the same configurations are showed in table (5):
α2 1 2 3 4 5 -0.17 -0.17 -0.17 -0.17 -0.17
α3 -0.28 -0.28 -0.28 -0.28 -0.28
Fitness 7.19E-06 7.19E-06 7.19E-06 7.19E-06 7.19E-06 6 7 8 9 10
α2 -0.47 -0.03 -0.17 -0.17 -0.17
α3 0.03 -0.38 -0.28 -0.28 -0.28
Fitness 3.59E-05 4.49E-05 7.19E-06 7.19E-06 7.19E-06
Table 5: Results of the first study Table (6) shows the comparison between the best results founded (smallest fitness or error) and the values of the original model.
k2 [N/m] 1.25E+09 Initial values Values after the update 1.04E+09 k3 [N/m] 1.43E+09 1.03E+09 First natural frequency [Hz] 5.14 4.57 Fitness 0.12 7.19E-06
Table 6: Comparison between the initial values and the best results in the 1st study The results illustrate that 8 simulations with the genetic algorithms found the same values, k2 was set at -17% and k3 by -28%, with an error of 7,19x10-6. Figure (16) helps in assessing the efficiency of the algorithm by viewing the surface of the objective function of optimization problem with the position of global minimum.
0.5 0.4 0.3 0.25 0.2 0.1 0.2 0.3
α3
0 0.15 ?0.1 ?0.2 ?0.3 0.05 ?0.4 ?0.5 ?0.5 Global minimum
X= ?0.17 Y= ?0.28 Level= 7.1867e?006
0.1
0
0.5
α
2
Figure 16: Objective function of the problem in the first study indicating the global minima. After analyzing table (6) and fig. (16) it can be concluded that the algorithm for optimization achieved the global minimum (α2 = -0.17; α3 = -0.28) in 80% of the cases, despite the large range with many local minima (dark blue region). This result shows the good ability of heuristic algorithms to find optimal solutions to complex problems. 7.4.2 Second study: multiobjective optimization with simulated vibration data
In the second study were used the natural frequencies identified with the data from integration of the motion equation of the hydrogenerator (see section 7.2). In this case, one has a multiobjective optimization because the objective function consider simultaneously two natural frequencies. The parameters to be corrected are k2 and k5 that are the most sensitive to changes of the first two natural frequencies. The objective function used in this study is based on the weighted sum of the errors of the two frequencies using the methodology described in Section 6.1. This type of evaluation function was used because of the facility to implement and analyze the results. The objective function can be written as:
Min J
( p (α ) ) = w1 ×
ω1adj ( p(α ) ) ? ω1 exp ω1 exp
+ w2 ×
ω 2 adj ( p(α ) ) ? ω 2 exp , for -50% < α < 50% ω 2 exp
(16)
where ω1adj and ω2adj are the natural frequencies obtained by the solution of a problem of eigenvalue and eigenvector with the modified stiffness values, ω1exp and ω2exp are the natural frequencies identified with simulated vibration data, w1 and w2 is the weight coefficients, p is the vector with the k2 and k5 stiffness values that will be updated, α is the vector with the α2 and α5 correction values and J is the objective function. The algorithm genetic used to this problem used the same configurations of the Section 7.4.1 showed in Table 4 and the weights given to the both frequency errors used was equal: w1=w2=0.5. For this study, it was simulated a loss of stiffness of 5 and 10 % in k2 and k5 to observe if the genetic algorithm could restore the original values of stiffness used to simulate the “experimental” data. The results of 10 simulations for the two cases are shown in the table (7).
α2 0.055 0.055 0.055 0.055 0.055 α2 0.114 0.114 0.113 0.114 0.114 α5 0.053 0.053 0.053 0.053 0.053 α5 0.112 0.112 0.112 0.112 0.112 Loss of 5 % in k2 and k5 Fitness α2 4.88E-05 6 0.055 0.055 4.88E-05 7 4.88E-05 8 0.055 4.88E-05 9 0.055 4.88E-05 10 0.055 Loss of 10 % in k2 and k5 Fitness α2 6.17E-05 6 0.114 6.17E-05 7 0.114 9.31E-05 0.113 8 6.17E-05 9 0.114 6.17E-05 10 0.114 α5 0.053 0.054 0.053 0.053 0.053 α5 0.112 0.112 0.112 0.112 0.112 Fitness 4.88E-05 1.08E-04 4.88E-05 4.88E-05 4.88E-05 Fitness 6.17E-05 6.17E-05 9.31E-05 6.17E-05 6.17E-05
1 2 3 4 5
1 2 3 4 5
Table 7: Results of simulations with genetic algorithms for losses of 5 and 10% of stiffness Table (8) illustrate a comparison between the initial values of stiffness parameters, natural frequencies and the fitness of this configuration and of the best updated values of the 10 simulations.
k2 [N/m] 1.19E+09 1.25E+09 k2 [N/m] 1.13E+09 1.25E+09 Loss of 5% in k2 and k5 k5 [N/m] First and second natural frequencies 3.04E+10 [5.08 6.53] 3.20E+10 [5.14 6.63] Loss of 10% in k2 and k5 k5 [N/m] First and second natural frequencies 2.88E+10 [5.02 6.42] 3.20E+10 [5.14 6.63] Fitness 0.014 4.88E-05 Fitness 0.028 6.17E-05
Initial values Values after updating
Initial values Values after updating
Table 8: Comparison between the initial values and the best values founded in the 2nd study Table (7) and (8) shows that the optimization algorithm achieved the goal to recover the original values of k2 and k5 described in Section 2 and was also able to set the natural frequencies identified with the simulated data. Moreover, it was found a good level of repeatability in the results. To verify if the genetic algorithm reached the global minimum point, it was performed a process of exhaustive search where all the possible combinations to the problem are calculated and compared with the responses found by the genetic algorithm. Thus, it was concluded that the best values founded in the cases of loss of 5 and 10% of stiffness are the global minimum, as illustrate in fig (17).
Loss of 5% of stiffness
0.5 0.16 0.4 0.3 0.2 0.1
5
Loss of 10% of stiffness
0.5 0.16 0.4 0.14 0.12 0.1 0.3 0.2 0.1 0.14 0.12 0.1 0.08 0.06 ?0.2 0.04 0.02 ?0.3 ?0.4 ?0.5 ?0.5 0 0.5 0.04 0.02
α
0 ?0.1
α5
0.08 0.06
0 ?0.1
?0.2 ?0.3 ?0.4 ?0.5 ?0.5 0 0.5
Global minimum
α2
Global minimum
α2
a) b) Figure 17: Objective functions with the original stiffness with: a) 5% of loss and b) 10% of loss.
8
Final remarks
In this work it was proposed a model updating of a simplified model of a generating unit of Itaipu. In a first study it was used real vibration data to identify the first natural frequency contained in this signal ,with the ARMA model, and this value was used to update the model of the hidrogenerator with the objective to minimize the error between the experimental frequency and the frequency calculated with the updated model using genetic algorithms. In the second study, the analytic model was used to generate simulated vibration data of the radial and angular displacement and the two natural frequencies were identified as like in the first study. After this, the original model used to generate the simulated data was modified with a loss of stiffness and then, the updating procedure was made to verify if the genetic algorithm can retrieve the original values of the stiffness parameters in a multiobjective optimization problem. The modal analysis using the ARMA model presented real negative poles or complex conjugates where this is an index of stability of the model. The application of genetic algorithms in the process of model updating was successful in finding the global optimum of the problems validating the use of heuristic algorithms to inverse problems like model updating.
Acknowledgments
The second and the third author would like to thank the PTI - C&T and CNPq for their scholarships, respectively. The authors are thankful to the support provided by Parque Tecnológico Itaipu (PTI) and Itaipu Binacional.
References
[1] G. C. Brito Jr, “Application of Simplified Statistics Models in Hidro Generating Unit Health Monitoring”, “Damage Prognosis for Aerospace, Civil and Mechanical Systems”, John Wiley & Sons Ltd., West Sussex, England (2005). [2] J. E. Mottershead, M. I. Friswell, Model Updating in Structural Dynamics: a Survey, Journal of Sound and Vibration, Vol. 167, No. 2, (1993), pp. 347-375. [3] J. A. Pereira, W. Heylen, S. Lammens, P. Sas, Influence of the number of frequency point and resonance frequency of model Updating Techniques for Health Condition Monitoring and Damage
Detection in Flexible Structures, XII- International Modal Analysis Conference.- IMAC, Nashville, USA (1995). [4] G. Kim, Y. Park, An Improved Updating Parameter Selection Method and Finite Element Model Update using Multiobjective Optimization Technique, Mechanical Systems and Signal Processing , No. 18 (2004), pp. 59-78. [5] S. Zivanovic, A. Pavic, P. Reynolds, Finite Element Modeling and Updating of a Lively Footbridge: The Complete Process, Journal of Sound and Vibration, No. 301 (2007), pp. 126-145. [6] G. C. Brito Jr, K. A. Brol, K. B. Brol, Dynamic Behavior of Large Rotating Eletric Machines (Analysis with Simplified Analytical Models), International Symposium on Stability Control of Rotating Machinery, Calgary, Alberta, Canada (2007). [7] S. Silva,G. C. Brito Jr, Damage Detection in Large Rotating Electric Machines Using Time Series Analysis , 8? Brazilian Conference on Dynamics, Control and Applications (2009). [8] P. H. Wirsching, T. L. Paez and O. Heith, Random Vibrations: Theory and Practice, John Wiley & Sons, Inc (1995). [9] P. S. R. Diniz, E. A. B. Silva, S. Lima Netto, Digital Signal Processing. Bookman Company, 1st Edition (2004). [10] L. Ljung, System Identification: Theory for the Use, 2nd edition, Prentice-Hall (1998). [11] L. A. Aguirre, Introdu??o à Identifica??o de Sistemas - Técnicas Lineares e N?o-Lineares Aplicadas a Sistemas Reais. Editora UFMG, 3rd Edition (2007). [12] V. Steffen Jr, D. A. Rade, Model-Based Inverse Problems in Structural Dynamics, “Damage Prognosis for Aerospace, Civil and Mechanical Systems”, John Wiley & Sons Ltd., West Sussex, England (2005). [13] M. I. Friswell, J. E. Mottershead, Finite Element Model Updating in Structural Dynamics, Kluwer Academic Publishers, Dorbrecht , Netherlands (1995). [14] J. He, Z. Fu, Moda Analysis, Butterworth Heinemann, London, England (2001). [15] P. G. Bakir,E. Reynders, G. De Roeck, Sensitivity-based Finite Element Model Updating using Constrained Optimization with a Trust Region Algorithm, Journal of Sound and Vibration, No. 305 (2007), pp. 211-225. [16] M. C. Goldbarg, P. C. L. Luna, Otimiza??o Combinatória e Programa??o Linear – Modelos e Algoritmos, 1st Edition, Editora Campus, Rio de Janeiro, Brazil (2000). [17] A. Tebaldi, L. S. Coelho, V. Lopes Jr, Detec??o de falhas em estruturas inteligentes usando otimiza??o por nuvem de partículas: fundamentos e estudo de casos, SBA Controle e Automa??o, Vol. 17, No. 3 (2006), pp. 312-330 [18] R. I. Levin, N. A. J. Lieven, Dynamic Finite Element Model Updating using Simulated Annealing and Genetic Algorithms, Mechanical Systems and Signal Processing, Vol. 12 No. 1 (1997), pp. 91120. [19] D. E. Goldberg, “Genetic Algorithms in Search, Optimization & Machine Learning”, AddisonWesley, Massachusetts, United States (1989). [20] Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs, 3rd Edition, Springer, Charlotte, United States (1996). [21] J. B. Bodeux and J.C. Golinval, Application of ARMA models to the identi?cation and damage detection of mechanical and civil engineering structures. Smart Materials and Structures, Vol. 10, pp. 479-489 (2001). [22] J. Juang, M. Q. Phan, Identification and Control of Mechanical Systems, 1st edition, Cambridge University Press (2006).