# HOC Based Blind Identification of Hydroturbine Shaft Volterra System

Authors: Bing Bai; Lixiang Zhang
Journal: Shock and Vibration
 (2017)

---

## Abstract

In order to identify the quadratic Volterra system simplified from the hydroturbine shaft system, a blind identification method based on the third-order cumulants and a reversely recursive method are proposed. The input sequence of the system under consideration is an unobservable independent identically distributed (i.i.d.), zero-mean and non-Gaussian stationary signal, and the observed signals are the superposition of the system output signal and Gaussian noise. To calculate the third-order moment of the output signal, a computer loop judgment method is put forward to determine the coefficient. When using optimization method to identify the time domain kernels, we combined the traditional optimization algorithm (direct search method) with genetic algorithm (GA) and constituted the hybrid genetic algorithm (HGA). Finally, according to the prototype observation signal and the time domain kernel parameters obtained from identification, the input signal of the system can be gained recursively. To test the proposed method, three numerical experiments and engineering application have been carried out. The results show that the method is applicable to the blind identification of the hydroturbine shaft system and has strong universality; the input signal obtained by the reversely recursive method can be approximately taken as the random excitation acted on the runner of the hydroturbine shaft system.


---

## Body

## 1. Introduction

There are lots of nonlinear phenomena in engineering. Knowing of the coherent characteristics is of great importance to engineering safety. Thus exploring an appropriate method to describe an input-output relation defined by the nonlinear system has been a hot research topic. At present, Volterra series is a commonly used technique in researchers, and achievements are involved with plenty of engineering fields. Maheswaran and Khosa [1, 2] use a Volterra series model with wavelet functions to identify the nonlinearity of a natural stream flow and to forecast the water levels in groundwater system, respectively. In forensic investigation, a filter designed with Volterra series can well work when it is necessary to identify the latent fingerprints [3]. A nonlinear phenomenon that is produced by a distorted harmonic audible signal is recognized by estimating sound pressure level with Volterra series [4]. In addition, the nonlinear structure response of a large container carrier in irregular seaway is calculated with a third-order Volterra model [5]. For a multiple input-multiple output (MIMO) system, a topological assemblage scheme is developed to make synthesis of Volterra system and applied to engineering, such as vibration analysis of wind-excited suspended cable [6].The applications of Volterra series are mainly demonstrated in two aspects, that is, modeling and identification of an input-output system. Two important properties of Volterra series are (1) output of a Volterra system depends linearly on kernel parameters of this system; (2) nonlinearity of a signal can be represented throughout multidimensional operators on product of samples. Therefore, Volterra model is generally used to describe the input-output relationship of a nonlinear system in a long term. Nonlinear systems, such as system with additive white Gaussian noises [7], structure subjected to impact load [8], and turbomachinery [9], can be represented by Volterra series. Systems with multispectral properties are also reported [10, 11].Peng et al. make a lot of achievements on the modeling, identification, and prediction with Volterra series [12–14]. They discussed lots of problems, for instance, the influence on the force transmissibility of MDOF structures with cubic nonlinear damping device [12], a new method to detect the breathing fatigue crack in nonlinear beam [13] and the analytical expression of calculating the bispectrum with higher order cumulants for nonlinear systems subjected to zero-mean Gaussian excitation [14]. These works highly prove the advantages of Volterra series in describing a nonlinear system and the convenience for derivation of the subsequent formulae.It is fundamental to use Volterra series method to identify a nonlinear system. Take Duffing oscillator as an example, its “jump” phenomenon can be recognized by the Volterra series method [15], which is expanded to identify a bilinear oscillator [16]. To deal with identification problems more effectively, Volterra series is usually combined with an optimization method [17]. However, the input signal of a nonlinear system is not always known. Thus the blind identification method without the input signal is developed. Le Caillec [18] blindly identified a second-order Wiener model driven by a Gaussian noise with corrupted output data. Ghandchi Tehrani et al. [19] compared the Volterra series method with the harmonic balance method, which proves the ability of blind identification by the Volterra series method.Actually, when Volterra series is employed to model and blindly identify a nonlinear system, it needs to often meet the situation in which higher order cumulants are used. Antari et al. [20, 21] considered the third-order and fourth-order cumulants to blindly identify a Hammerstein system. However, the expression of higher cumulants sometimes seems to be quite complicated [22]. For a hydroturbine main shaft system in a hydropower station, the nonlinear input-output relation is very complicated. Especially for the input signals, such as the hydraulic excitation on the runner, its working mechanisms are not completely known so far. On the contrary, the corresponding output signals, such as the throw at turbine guide bearing, are always easy to be obtained. Thus it is a mathematical problem which needs to be solved frequently in hydraulic engineering to identify the whole shaft system with the observed vibration response only.This paper proposes a novel blind identification method based on higher order cumulants for a nonlinear system described by a hydroturbine shaft Volterra model. A computer loop judgment method is used to determine the coefficientsγijk when the third-order moment of the input signal is calculated. It greatly simplifies formula and algorithm. Finally, the Volterra kernel parameters of the nonlinear system under consideration are identified by the hybrid genetic algorithm (HGA).


## 2. Description of System

The shaft system of a certain hydroturbine generating set can be simplified as a quadratic Volterra model, as shown in Figure1. The input sequence of the system is xt, representing the hydraulic excitation on the runner; the output sequence yt denotes the dynamic displacements at turbine guide bearing; the observed signal is st, and Gaussian noise is gt.Figure 1
Hydroturbine shaft system and quadratic Volterra system.The corresponding nonlinear system is described as(1)yt=∑i=0q∑j=0qhi,jxt-ixt-j,(2)st=yt+gt,where t stands for the discrete time, hi,j stands for the time domain kernels of the quadratic nonlinear system, and q represents the system memory length q+1, and i, j are the time delay. The excitation sequence xt of the system is an unobservable, independent identically distributed (i.i.d.), zero-mean and non-Gaussian stationary signal. It is assumed that the following equations hold:(3)γ2x=Ex2t≠0,γ3x=Ex3t≠0,γ4x=Ex4t≠0,γ6x=Ex6t≠0,where γm,x stands for the mth-order moment at the origin of the input signal. It is proved that, when xt is a zero-mean sequence, the following properties [22] hold:(4)Ext1⋯xt6=γ6xift1=t2=t3=t4=t5=t6γ4xγ2xift1=t2≠t3=t4=t5=t6γ3x2ift1=t2=t3≠t4=t5=t6γ2x3ift1=t2≠t3=t4≠t5=t60otherwise.LettingCm,yτ1,…,τm-1 and Cm,sτ1,…,τm-1 be cumulants of the output sequence and the observed sequence, respectively, a relationship between the higher order cumulants (HOC) and time domain kernels hi,j of the quadratic nonlinear system is established by calculation. In the present study, we try to establish a relationship between the third-order cumulants of the observed signal and time domain kernels hi,j. Because the third-order or higher order cumulants of the Gaussian noise are zero, we obtain from (2) the following:(5)C3sτ1,τ2=C3yτ1,τ2+C3gτ1,τ2=C3yτ1,τ2.Therefore, the third-order cumulants of the output sequence yt need to only be calculated. For a stationary random sequence yt, the third-order cumulants are defined as(6)C3yτ1,τ2=M3yτ1,τ2-M1yM2yτ1+M2yτ2+M2yτ2-τ1+2M1y3,where M1y, M2yτ, and M3yτ1,τ2 denote the first-order moment, second-order moment, and third-order moment of the sequence yt, respectively. The specific expression of M1y is stated as(7)M1y=∑i=0q∑j=0qhi,jExt-ixt-j=γ2x·∑i=0qhi,i.It implies that when i,j are not equal, the expectation value is zero; only the diagonal elements hi,i are left for the time domain kernels hi,j. M2yτ is as follows:(8)M2yτ=Eytyt+τ=∑i1,i2,j1,j2=0qγij·hi1,i2hj1,j2,where γij is the coefficient, which consists of the second-order moment or the fourth-order moment of the input signal xt, and the value depends on the different time delay of the sequence xt. Among the four subscripts -i1,-i2,τ-j1, and τ-j2, if the value of one subscript (time) is different from the rest of the three, the expectation is zero, and γij=0. The values of γij are yielded as(9)γij=γ4xif-i1=-i2=τ-j1=τ-j2γ2x2if-i1=-i2≠τ-j1=τ-j2,orif-i1=τ-j1≠-i2=τ-j2,orif-i1=τ-j2≠-i2=τ-j10otherwise.In order to calculate the three second-order cumulants in (6), we substitute τ, respectively, by τ1, τ2, and τ2-τ1 in (8). The expression of the last term M3yτ1,τ2 in (6) is written as(10)M3yτ1,τ2=Eytyt+τ1yt+τ2=∑i1,i2,j1,j2,k1,k2=0qγijk·hi1,i2hj1,j2hk1,k2in which coefficient γijk is composed of one or several certain order moments of the input signal sequence xt. According to the different time delay of the sequence xt, the different values are obtained as (11)γijk=γ6xif-i1=-i2=τ1-j1=τ1-j2=τ2-k1=τ2-k2γ4xγ2xif-i1=-i2≠τ1-j1=τ1-j2=τ2-k1=τ2-k2γ3x2if-i1=-i2=τ1-j1≠τ1-j2=τ2-k1=τ2-k2γ2x3if-i1=-i2≠τ1-j1=τ1-j2≠τ2-k1=τ2-k20otherwise.The judgment of the first line γ6x is simple, and the judgments of the rest γ4xγ2x, γ3x2, and γ2x3 are relatively complicated since there exist a variety of different combinations, and all combinations need to compile a programming code to loop for judgment, which is called as a method on the computer loop judgment (details in the Appendix). With this method, it is possible to calculate the second-order moment and third-order moment of the output sequence yt and to have a great superiority in computing.


## 3. Hybrid Genetic Algorithm

### 3.1. Consistent Estimate

For a known observed signalst, we need to have an expectation operation. Because the number of the observed signals is limited, it is difficult to precisely calculate the expectation required. So, an approximating method is used to calculate this value by consistent estimate; namely, (12)C^3sτ1,τ2=M^3sτ1,τ2-M^1sM^2sτ1+M^2sτ2+M^2sτ2-τ1+2M^1s3,(13)M^3sτ1,τ2=1N∑n=1Nsnsn+τ1sn+τ2,(14)M^2sτ=R^sτ=1N∑n=1Nsnsn+τ,(15)M^1s=μ^s=1N∑n=1Nsn.In order to complete the consistent estimate of the third-order cumulants in (12), we replace τ with τ1, τ2, and τ2-τ1 in (14), respectively.


### 3.2. Algorithm

The hybrid genetic algorithm (HGA) is an intelligent optimization algorithm which combines the simple genetic algorithm (SGA) with the traditional optimization methods.As we know, for a kind of optimization algorithm, the global search ability of the genetic algorithm (GA) is very strong. However, when the search process is close to a local optimal point, its local search ability starts to become poor. That is a shortcoming of the simple genetic algorithm. On the contrary, the traditional direct search method of discrete variables is lack of global search ability, although the local search ability is very strong. Therefore, the authors make a hybrid genetic algorithm by combining the traditional direct search method with the genetic algorithm. The hybrid genetic algorithm (HGA) improves the ability of the global and local searches in calculation.The basic steps of HGA are sentenced in the process of GA. The local search is carried out in a certain probabilityph (for example 5%) to find out a local optimal point as soon as possible. Then it goes back to the genetic operation for the global search to complete the global search. The flow chart of HGA is as shown in Figure 2.Figure 2
Flow chart of hybrid GA.


### 3.3. Blind Identification of Time Domain Kernels

Through the above derivation, we calculate the consistent estimates of the third-order cumulants of the observed signalst by using (12) and (13)–(15). On the other hand, we obtain a relationship between time-domain kernels hi,j of the nonlinear system and the third-order cumulants by using (6) and (7)–(11). Thus, letting (6) equal (12), we obtain an equation as(16)C^3sτ1,τ2=M3yτ1,τ2-M1yM2yτ1+M2yτ2+M2yτ2-τ1+2M1y3.Obviously, this is a cubic nonlinear equation with unknown hi,j. By solving this equation, hi,j is obtained. It is noted that because the output signal of the system is only used in calculation, namely, no input signal xt is needed, so, this identification method is called blind identification.In actual application, the third-order cumulants ofq+12 have to be computed; that is, C3yτ1,τ2, 0≤τ1, τ2≤q. Since it is not easy to directly solve a cubic nonlinear system, the above-mentioned HGA is proposed to solve (16). The objective function (the fitness value of an individual) of the optimization problem is as follows:(17)min⁡Fh=∑τ1=0q∑τ2=τ1qC3yτ1,τ2-C^3sτ1,τ22,where h is the unknown column matrix composed of the elements in time domain kernels hi,j. Due to the symmetry properties of the third-order cumulants, namely, C3yτ1,τ2=C3yτ2,τ1, the loop variable τ2 in (17) increases from τ1 to q.For HGA, when the Gaussian noise is zero-mean i.i.d. signal, the following constraint condition can be used:(18)M1y-M^1s≤ε,where ε denotes the admissible error of the constraint condition. Comparing (7) with (15), it is seen that the introduction of the constraint (18) can guarantee the small deviation of the diagonal elements hi,i in kernels hi,j to be identified.


## 3.1. Consistent Estimate

For a known observed signalst, we need to have an expectation operation. Because the number of the observed signals is limited, it is difficult to precisely calculate the expectation required. So, an approximating method is used to calculate this value by consistent estimate; namely, (12)C^3sτ1,τ2=M^3sτ1,τ2-M^1sM^2sτ1+M^2sτ2+M^2sτ2-τ1+2M^1s3,(13)M^3sτ1,τ2=1N∑n=1Nsnsn+τ1sn+τ2,(14)M^2sτ=R^sτ=1N∑n=1Nsnsn+τ,(15)M^1s=μ^s=1N∑n=1Nsn.In order to complete the consistent estimate of the third-order cumulants in (12), we replace τ with τ1, τ2, and τ2-τ1 in (14), respectively.


## 3.2. Algorithm

The hybrid genetic algorithm (HGA) is an intelligent optimization algorithm which combines the simple genetic algorithm (SGA) with the traditional optimization methods.As we know, for a kind of optimization algorithm, the global search ability of the genetic algorithm (GA) is very strong. However, when the search process is close to a local optimal point, its local search ability starts to become poor. That is a shortcoming of the simple genetic algorithm. On the contrary, the traditional direct search method of discrete variables is lack of global search ability, although the local search ability is very strong. Therefore, the authors make a hybrid genetic algorithm by combining the traditional direct search method with the genetic algorithm. The hybrid genetic algorithm (HGA) improves the ability of the global and local searches in calculation.The basic steps of HGA are sentenced in the process of GA. The local search is carried out in a certain probabilityph (for example 5%) to find out a local optimal point as soon as possible. Then it goes back to the genetic operation for the global search to complete the global search. The flow chart of HGA is as shown in Figure 2.Figure 2
Flow chart of hybrid GA.


## 3.3. Blind Identification of Time Domain Kernels

Through the above derivation, we calculate the consistent estimates of the third-order cumulants of the observed signalst by using (12) and (13)–(15). On the other hand, we obtain a relationship between time-domain kernels hi,j of the nonlinear system and the third-order cumulants by using (6) and (7)–(11). Thus, letting (6) equal (12), we obtain an equation as(16)C^3sτ1,τ2=M3yτ1,τ2-M1yM2yτ1+M2yτ2+M2yτ2-τ1+2M1y3.Obviously, this is a cubic nonlinear equation with unknown hi,j. By solving this equation, hi,j is obtained. It is noted that because the output signal of the system is only used in calculation, namely, no input signal xt is needed, so, this identification method is called blind identification.In actual application, the third-order cumulants ofq+12 have to be computed; that is, C3yτ1,τ2, 0≤τ1, τ2≤q. Since it is not easy to directly solve a cubic nonlinear system, the above-mentioned HGA is proposed to solve (16). The objective function (the fitness value of an individual) of the optimization problem is as follows:(17)min⁡Fh=∑τ1=0q∑τ2=τ1qC3yτ1,τ2-C^3sτ1,τ22,where h is the unknown column matrix composed of the elements in time domain kernels hi,j. Due to the symmetry properties of the third-order cumulants, namely, C3yτ1,τ2=C3yτ2,τ1, the loop variable τ2 in (17) increases from τ1 to q.For HGA, when the Gaussian noise is zero-mean i.i.d. signal, the following constraint condition can be used:(18)M1y-M^1s≤ε,where ε denotes the admissible error of the constraint condition. Comparing (7) with (15), it is seen that the introduction of the constraint (18) can guarantee the small deviation of the diagonal elements hi,i in kernels hi,j to be identified.


## 4. Numerical Experiment

To examine the proposed method, two quadratic nonlinear Volterra systems, E1 by (19) and E2 by (20), are herein considered as(19)yn=x2n-0.30x2n-1-0.85x2n-2+0.66x3n-3,(20)yn=x2n-0.3xnxn-1+0.5xnxn-2-0.26xnxn-3-0.30x2n-1+0.36xn-1xn-2-0.28xn-1xn-3-0.85x2n-2-0.2xn-2xn-3+0.66x2n-3.In (19) and (20), n is used to replace time t, to denote discrete time.Because (19) is of completely squared terms, the coefficients of each term are the true solutions. Therefore, the true solutions of E1 are(21)hi,i=1.00-0.30-0.850.66.In (20), the coefficient of a squared term is the true solution; however, half of the coefficient of a cross term is the true solution. Hence the true values of E2 are(22)hi,j=1.00-0.150.25-0.13-0.300.18-0.14-0.85-0.10sym0.66.The input signal x0n is generated by RAND command in MATLAB. Their lengths are 1024 and 8192, and each has 50 groups. Then the generated uniform white noises are scaled as independent identically distributed (i.i.d.) signal with variance 1. Assuming that the memory length of the model is M=q+1=4, the simulation is under the condition of 50 Monte-Carlo runs. At the same time, the different values of the signal-to-noise ratio (SNR) defined by (23) are considered [20]:(23)SNR=10·log10Ey2tEg2t.Only the diagonal elements exist in the time domain kernels hi,i of E1; in other words, the system is a simplified Volterra Model: Hammerstein Model. The identification results are listed in Table 1 for four cases of SNR=0,20,40,∞. At the same time, to assess the universality of the method, the time domain hi,i is identified by treating it as a full matrix with element hi,j. The results are listed in Table 2. The time domain kernels hi,j of E2 are originally full matrix with the elements, and the E2’s results are listed in Table 3.Table 1
Identified results [E1: diagonal elementhi,i].
True values SNR = 0 dB SNR = 20 dB SNR = 40 dB SNR =∞ Mean± SD Mean± SD Mean± SD Mean± SD h 0,0 1.00 0.9437 ± 0.1250 0.9877 ± 0.0913 0.9930 ± 0.0857 0.9962 ± 0.0685 h 1,1 −0.30 - 0.2460 ± 0.1368 - 0.3018 ± 0.1124 - 0.2997 ± 0.1138 - 0.3010 ± 0.1145 h 2,2 −0.85 - 0.7958 ± 0.1394 - 0.8328 ± 0.1003 - 0.8377 ± 0.1001 - 0.8390 ± 0.0872 h 3,3 0.66 0.6053 ± 0.1356 0.6557 ± 0.0981 0.6540 ± 0.0990 0.6534 ± 0.0953 F  0.751054 E - 01 0.316731 E - 02 0.314600 E - 02 0.357726 E - 02Table 2
Identified results [E1: full elementhi,j].
True values SNR = 0 dB SNR = 20 dB SNR = 40 dB SNR =∞ Mean± SD Mean± SD Mean± SD Mean± SD h 0,0 1.00 0.9568 ± 0.0631 0.9745 ± 0.0575 0.9783 ± 0.0610 0.9731 ± 0.0501 h 0,1 0 - 0.0161 ± 0.0768 - 0.0162 ± 0.0712 0.0001 ± 0.0712 - 0.0071 ± 0.0717 h 0,2 0 0.0054 ± 0.0658 - 0.0038 ± 0.0602 - 0.0045 ± 0.0645 - 0.0083 ± 0.0624 h 0,3 0 - 0.0147 ± 0.0426 - 0.0061 ± 0.0501 - 0.0162 ± 0.0517 - 0.0186 ± 0.0533 h 1,1 −0.30 - 0.2752 ± 0.0687 - 0.2928 ± 0.0706 - 0.2935 ± 0.0733 - 0.2948 ± 0.0668 h 1,2 0 - 0.0484 ± 0.0785 - 0.0212 ± 0.0693 - 0.0405 ± 0.0600 - 0.0294 ± 0.0635 h 1,3 0 0.0012 ± 0.0676 - 0.0012 ± 0.0579 - 0.0063 ± 0.0597 - 0.0047 ± 0.0570 h 2,2 −0.85 - 0.8289 ± 0.0668 - 0.8324 ± 0.0678 - 0.8348 ± 0.0713 - 0.8328 ± 0.0598 h 2,3 0 - 0.0096 ± 0.0782 - 0.0191 ± 0.0627 - 0.0143 ± 0.0579 - 0.0044 ± 0.0632 h 3,3 0.66 0.6448 ± 0.0604 0.6543 ± 0.0627 0.6528 ± 0.0639 0.6534 ± 0.0516 F  0.693044 E - 01 0.201256 E - 02 0.184063 E - 02 0.259258 E - 02Table 3
Identified results [E2: full elementhi,j].
True values SNR = 0 dB SNR = 20 dB SNR = 40 dB SNR =∞ Mean± SD Mean± SD Mean± SD Mean± SD h 0,0 1.00 0.9793 ± 0.0689 1.0118 ± 0.0523 1.0138 ± 0.0546 1.0191 ± 0.0507 h 0,1 −0.15 - 0.1713 ± 0.0945 - 0.1395 ± 0.0593 - 0.1364 ± 0.0543 - 0.1319 ± 0.0518 h 0,2 0.25 0.2185 ± 0.1070 0.2273 ± 0.0789 0.2324 ± 0.0768 0.2282 ± 0.0773 h 0,3 −0.13 - 0.0719 ± 0.0953 - 0.0890 ± 0.0604 - 0.0907 ± 0.0609 - 0.0972 ± 0.0564 h 1,1 −0.30 - 0.3062 ± 0.0699 - 0.3024 ± 0.0604 - 0.3071 ± 0.0620 - 0.3167 ± 0.0599 h 1,2 0.18 0.2247 ± 0.0913 0.1723 ± 0.0654 0.1765 ± 0.0580 0.1812 ± 0.0529 h 1,3 −0.14 - 0.0892 ± 0.1079 - 0.1223 ± 0.0693 - 0.1246 ± 0.0704 - 0.1183 ± 0.0685 h 2,2 −0.85 - 0.8162 ± 0.0606 - 0.8569 ± 0.0660 - 0.8494 ± 0.0662 - 0.8467 ± 0.0603 h 2,3 −0.10 - 0.0878 ± 0.0777 - 0.1008 ± 0.0520 - 0.1014 ± 0.0518 - 0.1001 ± 0.0513 h 3,3 0.66 0.6423 ± 0.0722 0.6586 ± 0.0707 0.6558 ± 0.0635 0.6567 ± 0.0611 F  0.107841 E + 00 0.351261 E - 02 0.340253 E - 02 0.384778 E - 02It is seen from Tables1–3 that the proposed method is suitable to well identify the quadratic nonlinear system under consideration. Table 1 shows that, when SNR = 0 dB, the relative errors of the four values of the time domain kernels hi,i are basically within 10%, and the standard deviation is between 0.12 and 0.14. With the increase of the signal-to-noise ratio, the relative errors are gradually reduced to less than 2%, and the standard deviation is to 0.06–0.11; the objective function value (sum of the square of the third-order cumulants errors) is also reduced. The identification error becomes large with the small SNR; this is due to the fact that there exists certain error between Gaussian noise signal generated by the numerical method and the real Gaussian distribution, which leads a direct effect of the identification results.In Table2, the same sequence to the signal data as Table 1 is used. The identified mean values of the nondiagonal elements are generally under 0.02 (the real value is 0), and the standard deviation is below 0.08. This is very important; when the blind identification method is used to identify a quadratic system, we do not know the length of the memory, or whether the time domain kernels hi,j have only the diagonal elements. By using different memory length to identify an identical system, we obtained very similar results, which show that this method has strong universality.


## 5. Blind Identification of Quadratic Volterra System of Number 4 Hydroturbine Shaft in a Certain Hydropower Station

The number 4 hydroturbine generating set shaft system in a certain hydropower station, as shown in Figure1(a), can be simplified as a second-order Volterra system (shown in Figure 1(b)). The observed data of the lateral throw y2n at the turbine guide bearing in the experiment is shown in Figure 3 (after filtration). These observed signals are taken as the output signals yn, and this engineering problem is identified according to the above proposed method. Note that, for the problem in which the input signals are unknown, there is a precondition here: the input must be zero-mean, i.i.d., stationary signals, and their statistics are known.Figure 3
Transient waveshape of the lateral throw at turbine guide bearing (filtered).According to the 50 groups of the zero-mean, i.i.d. stationary signals generated by MATLAB, the statistics of the input signals are chosen as follows:(24)γ2x=Ex2n=1.000,γ3x=Ex3n=2.645×10-2,γ4x=Ex4n=1.800,γ6x=Ex6n=3.900.The memory length of the system is taken as M2=q2+1=6, so the time domain kernels parameters (optimization variables) to be identified are(25)hi,j=h0,0h0,1h0,2h0,3h0,4h0,5h1,1h1,2h1,3h1,4h1,5h2,2h2,3h2,4h2,5h3,3h3,4h3,5h4,4h4,5h5,5.The system becomes an optimization problem with 21 design variables. When using HGA, we set the population as 50, crossover probability 50%, mutation probability 5%, and local search probability 5%. The lower and upper bounds of the variables are[-100,100], and the increment of the discrete variables is 0.01. The interval of each variable is narrowed down after some test run, and the increment of the discrete variables is ultimately decreased to 0.00001. The number of the evolution generations is one million.After identification, the time domain kernels parametershi,j of the Volterra system for this shaft system are obtained, as shown in Figure 4. The frequency domain kernels parameters (Generalized Frequency Response Function, GFRF) can be obtained by two-dimensional Fourier Transform, as shown in Figure 5.Figure 4
The second-order time domain kernels of the shaft system in overspeed test.Figure 5
The frequency domain kernels of the shaft system in overspeed test.So far, by use of the blind identification method proposed in this paper, the blind identification of the system was successfully carried out, based on the observed throw signals at the turbine guide bearing.


## 6. Reversely Recursive Solution to Hydraulic Excitation

### 6.1. Reversely Recursive Algorithm

The hydraulic excitation acted on the hydraulic turbine runner is very important in the design and calculation of shaft system in a hydroelectric generating set but difficult to be observed. A reversely recursive algorithm is proposed in the paper, which can be used to approximately find out the input signal of the system based on the following two sets of data: (1) the observed output signal y0n of the system; (2) the system time domain kernels hi,j obtained through blind identification.Expanding the second-order discretized time domain Volterra system (1), the specific expression of yn can be obtained. When n=0, from the following expression:(26)y0=h0,0x0x0.We can get(27)x0=±y00h0,0,where y00 is the observed output signal of the system. x0 thus obtained has two roots.Whenn=1,(28)y1=h0,0x1x1+h0,1x1x0+h1,0x0x1+h1,1x0x0in which x0 is known and x1 is unknown. Equation (28) becomes a quadratic equation of x1:(29)h0,0x1x1+h0,1x0+h1,0x0x1+h1,1x0x0-y1=0.Set a=h0,0, b=h0,1x0+h1,0x0=2h0,1x0, and c=h1,1x0x0-y01; then(30)a·x1x1+b·x1+c=0.So(31)x1=-b±b2-4ac2a.Becausex0 has two roots, substituting, respectively, the two roots into (31), we can obtain the four roots of x1. The question now is which roots should x0 and x1 take? The answer is obvious. First, calculate y0 and y1 through (26) and (28), then calculate the errors of y0-y00 and y1-y01, and finally compare the errors. The roots with minimum error are the true solution. So we take(32)F=y1-y012+y0-y002as the optimal objective function and find out the minimum value Fmin from several values of objective function F, and the corresponding x0 and x1 are taken as the true solution (input signal).When the previous values ofxn-2 and the x· have been confirmed, the general formulae for calculating xn are as follows:(33)yn=h0,0xn-0xn-0+h0,1xn-0xn-1+⋯+h0,qxn-0xn-q+h1,0xn-1xn-0+h1,1xn-1xn-1+⋯+h1,qxn-1xn-q+⋯+hq,0xn-qxn-0+hq,1xn-qxn-1+⋯+hq,qxn-qxn-q.Leta=h0,0, b=∑i=1q2·h0,ixn-i, and c=∑i=1q∑j=1qhi,jxn-ixn-j-y0n; thus(34)xn=-b±b2-4ac2a.Becausexn-1 has two values, we can obtain the four roots of xn. In order to determine xn-1 and xn, we still use optimization method. The objective function for optimization is(35)F=∑i=0qyn-i-y0n-i2.Finding the minimum valueFmin, the corresponding xn-1 and xn can be taken as the input signal. Thus through recursive calculation, we can find out x0~xN; here N+1 is the sampling length. The above method is the reversely recursive algorithm.


### 6.2. Numerical Experiment  E2

For convenience, we still use the experimental data E2 in Section4. According to the above method, the calculated input data xn can be obtained, as shown in Figure 6, where xx(n) is the calculated result and x0(n) is the original input signal.Figure 6
E2 numerical experiment results (input signals).To illustrate the correctness and feasibility of the calculation results, substituting in (1) the calculated results xx(n) and the original input signal x0(n), respectively, we can find out the output signal y(n) and the original output signal y0(n) and their differences:(36)dyn=yn-y0nand draw in Figure 7. At the same time, the relative error of output signal energy is calculated by using the following formula:(37)RE=Σy2n-Σy02nΣy02n×100%.Figure 7
The comparison of the calculated output signals to the original observed signals.From Figures6 and 7 we can see that, by the reversely recursive method proposed in this paper, we can approximately find out the input signal of the system based on the system output signal. The calculated input signal xx(n) is slightly different from the original input signal x0(n), which is mainly caused by the quadratic nonlinearity of the system. However, the errors dy(n) between the system output signal y(n) calculated through the input signal xx(n) and the original output signal y0(n) are smaller (Figure 7). Only a few signal values have errors, and most output signal values are the same as the original value. The relative error of the output signal energy by (37) is only 3.96%. Therefore, the input data xx(n) thus calculated can be approximately used as the system input signal.


### 6.3. Reversely Recursive Solution to Hydraulic Excitation on Turbine

The hydraulic excitation (input signal)xx(n) on the turbine runner and its distribution can be calculated recursively using the proposed method in Section 6.1, as shown in Figure 8.Figure 8
The input signals calculated by recursive method and their histogram.
(a)
Waveform (b)
HistogramSince there is no real input signals that can be compared, we cannot tell the difference between the input signalxx(n) obtained with recursive method and the real input signals. However, we can compare the calculated output signal yn with the real output signal y0n. Figure 9 shows the calculated output signal yn and the absolute errors with observed output signal y0n. From the figure we can see that the error is not very big. Besides, the relative error of the output signal energy was calculated, whose value is 8.95%. This shows that the method is feasible; that is to say, the input signal xx(n) is acceptable and can be approximately taken as the hydraulic excitation on the turbine runner of a hydroelectric generating set shaft system.Figure 9
The output signalsy(n) calculated by xx(n) and their errors.
(a)
y ( n ) waveform (b)
Errors compared with the observed signal


## 6.1. Reversely Recursive Algorithm

The hydraulic excitation acted on the hydraulic turbine runner is very important in the design and calculation of shaft system in a hydroelectric generating set but difficult to be observed. A reversely recursive algorithm is proposed in the paper, which can be used to approximately find out the input signal of the system based on the following two sets of data: (1) the observed output signal y0n of the system; (2) the system time domain kernels hi,j obtained through blind identification.Expanding the second-order discretized time domain Volterra system (1), the specific expression of yn can be obtained. When n=0, from the following expression:(26)y0=h0,0x0x0.We can get(27)x0=±y00h0,0,where y00 is the observed output signal of the system. x0 thus obtained has two roots.Whenn=1,(28)y1=h0,0x1x1+h0,1x1x0+h1,0x0x1+h1,1x0x0in which x0 is known and x1 is unknown. Equation (28) becomes a quadratic equation of x1:(29)h0,0x1x1+h0,1x0+h1,0x0x1+h1,1x0x0-y1=0.Set a=h0,0, b=h0,1x0+h1,0x0=2h0,1x0, and c=h1,1x0x0-y01; then(30)a·x1x1+b·x1+c=0.So(31)x1=-b±b2-4ac2a.Becausex0 has two roots, substituting, respectively, the two roots into (31), we can obtain the four roots of x1. The question now is which roots should x0 and x1 take? The answer is obvious. First, calculate y0 and y1 through (26) and (28), then calculate the errors of y0-y00 and y1-y01, and finally compare the errors. The roots with minimum error are the true solution. So we take(32)F=y1-y012+y0-y002as the optimal objective function and find out the minimum value Fmin from several values of objective function F, and the corresponding x0 and x1 are taken as the true solution (input signal).When the previous values ofxn-2 and the x· have been confirmed, the general formulae for calculating xn are as follows:(33)yn=h0,0xn-0xn-0+h0,1xn-0xn-1+⋯+h0,qxn-0xn-q+h1,0xn-1xn-0+h1,1xn-1xn-1+⋯+h1,qxn-1xn-q+⋯+hq,0xn-qxn-0+hq,1xn-qxn-1+⋯+hq,qxn-qxn-q.Leta=h0,0, b=∑i=1q2·h0,ixn-i, and c=∑i=1q∑j=1qhi,jxn-ixn-j-y0n; thus(34)xn=-b±b2-4ac2a.Becausexn-1 has two values, we can obtain the four roots of xn. In order to determine xn-1 and xn, we still use optimization method. The objective function for optimization is(35)F=∑i=0qyn-i-y0n-i2.Finding the minimum valueFmin, the corresponding xn-1 and xn can be taken as the input signal. Thus through recursive calculation, we can find out x0~xN; here N+1 is the sampling length. The above method is the reversely recursive algorithm.


## 6.2. Numerical Experiment  E2

For convenience, we still use the experimental data E2 in Section4. According to the above method, the calculated input data xn can be obtained, as shown in Figure 6, where xx(n) is the calculated result and x0(n) is the original input signal.Figure 6
E2 numerical experiment results (input signals).To illustrate the correctness and feasibility of the calculation results, substituting in (1) the calculated results xx(n) and the original input signal x0(n), respectively, we can find out the output signal y(n) and the original output signal y0(n) and their differences:(36)dyn=yn-y0nand draw in Figure 7. At the same time, the relative error of output signal energy is calculated by using the following formula:(37)RE=Σy2n-Σy02nΣy02n×100%.Figure 7
The comparison of the calculated output signals to the original observed signals.From Figures6 and 7 we can see that, by the reversely recursive method proposed in this paper, we can approximately find out the input signal of the system based on the system output signal. The calculated input signal xx(n) is slightly different from the original input signal x0(n), which is mainly caused by the quadratic nonlinearity of the system. However, the errors dy(n) between the system output signal y(n) calculated through the input signal xx(n) and the original output signal y0(n) are smaller (Figure 7). Only a few signal values have errors, and most output signal values are the same as the original value. The relative error of the output signal energy by (37) is only 3.96%. Therefore, the input data xx(n) thus calculated can be approximately used as the system input signal.


## 6.3. Reversely Recursive Solution to Hydraulic Excitation on Turbine

The hydraulic excitation (input signal)xx(n) on the turbine runner and its distribution can be calculated recursively using the proposed method in Section 6.1, as shown in Figure 8.Figure 8
The input signals calculated by recursive method and their histogram.
(a)
Waveform (b)
HistogramSince there is no real input signals that can be compared, we cannot tell the difference between the input signalxx(n) obtained with recursive method and the real input signals. However, we can compare the calculated output signal yn with the real output signal y0n. Figure 9 shows the calculated output signal yn and the absolute errors with observed output signal y0n. From the figure we can see that the error is not very big. Besides, the relative error of the output signal energy was calculated, whose value is 8.95%. This shows that the method is feasible; that is to say, the input signal xx(n) is acceptable and can be approximately taken as the hydraulic excitation on the turbine runner of a hydroelectric generating set shaft system.Figure 9
The output signalsy(n) calculated by xx(n) and their errors.
(a)
y ( n ) waveform (b)
Errors compared with the observed signal


## 7. Conclusions

A blind identification method based on the third-order cumulants is proposed for a quadratic Volterra nonlinear system. Through three numerical experiments and engineering application, the following conclusions are drawn.(1)
Two different types of the memory length are used to identify the same sets of the signal data sequence of hydroturbine main shaft, and the error between the two groups of the time domain kernelshi,j is small. It shows that the proposed method is applicable to the blind identification of the hydroturbine main shaft system and has strong universality.(2)
The proposed computer loop judgment method to determine the coefficientγijk for calculating the third-order moment of the output signal greatly simplifies the formulae derivation and compiling of the programming code and is suitable for calculating the fourth- or higher order cumulants.(3)
With the gradual increase of the signal-to-noise ratio, the relative error of the identified time domain kernels decreases, the standard deviation is reduced, and the objective function value is also reduced.(4)
According to the observed signal of the shaft system of a certain hydroturbine generating set, the time domain kernels of the Volterra system can be blindly identified by the method proposed in this paper. Furthermore, according to these time domain kernels and the observed signal, the input signal of the system can be recursively calculated. The results show that the computed input signal can be approximately taken as the stochastic excitation acted on the turbine runner, and the error is not big (relative error of the output signal energy is acceptable). The proposed method is of enlightening significance for the form of hydraulic action.



---

*Source: 6732704-2017-01-10.xml*