Qualitative study of Riccati difference equation on maneuvering target tracking and fault diagnosis of wind turbine gearbox

Abstract Various performance analyses using Interacting Multiple Model (IMM) algorithm for the stability of tracking a maneuvering target are dealt with the state error covariance. A noble way of analyzing the stability of the IMM algorithm for a linear systems with the upper and lower bounds of the error covariance is analyzed. For this persuasion, a two-model Interacting Multiple Model Kalman Filter (IMMKF) with constant acceleration and constant jerk model has been applied for two different case studies. One is tracking a maneuvering target, and the another is tracking a vibration of wind turbine gearbox, which helps to identify failure component in wind energy system. The required data are collected from a radar and a defected gear box of a test wind turbine, and the efficiency of IMM algorithm is analyzed by simulation experiments.


PUBLIC INTEREST STATEMENT
In this study, we have implemented the theoretical part of "Properties of Riccati difference equations" and provide the link between state error covariance of Kalman filter which is similar to Riccati difference equation. Since the Kalman filter has many technological applications and thus we have applied to two fields namely, target tracking in aircraft and wind turbine gearbox. In order to estimate the upper and lower bounds of state error covariance by using Riccati difference equation, we have used constant acceleration model and constant jerk model in the aircraft. These bounds are also used to analyze the stability of tracking.
In the wind turbine, the vibration of acceleration data was sensed by 12 sensors. Using Interacting Multiple Model Kalman Filter algorithm, upper and lower bounds of state error covariance were computed and this helps us to diagnose the fault in the gearbox.

Introduction
In numerous fields of engineering theory, especially in optimal control and filtering theory the Riccati equation enacts an essential role. Successful applications of this theory is witnessed in various practical engineering systems such as robot manipulators, underwater vehicles, power systems, electrical motors, automotive engines, etc. Nowadays, a contemporary way of solving the problems in discrete time non-linear systems using filter and optimal control has gained much attention. Generally, a filter is a collection of specified mathematical equations, which provide a dynamic computational means to determine the state of a process that reduces the mean square error. The support extended towards the estimation of current and previous states and also the forthcoming states is a powerful aspect of the filter. Estimation projection is an analytical method of incorporating state equality constraints in the Kalman filter (Berg, 1983). Here the filter is established just as the well-known relations among the state variables are satisfied by the state estimate. The determination of the steady state solution of the Riccati difference equation directly without iteration is discussed in (Vaughan, 1970) by D.R Vaughan and the robustness of fast computational for the discrete Riccati equation's solution is presented in (Lainiotis, 1975) by D.G. Lainiotis. Multiple filter model enables the accuracy of target state estimation with dynamics. The provision of a structure through the interacting multiple model (IMM) algorithm efficiently manages multiple filter models. The estimated states are generally the kinematic quantity of position, velocity, and acceleration (Bar-Shalom & Rong Li, 1995). Filters are utilized to reduce the uncertainty caused by noise in observation and to estimate quantities not directly observed in these measurements of condition. The filter uses a state process model that can be used to predict the behavior of the observed target accurately to estimate the desired kinematic quantities.
In the filtering problem, P 2 k À Á is physically interpreted as the prediction error covariance of Equation (2.12) in (Chan, Goodwin, & Sin, 2002). Consequently, is therefore a real, symmetric, and non-negative definite. A complete investigation of the problems and techniques of tracking targets is studied by Li and Jilkov in (Li & Bar-Shalom, 1991;Li & Jilkov, 2003 and a different approach of survey for selfadjusting IMM algorithm is contributed by Y.Bar-Shalom in (Mazor, Averbuch, Bar-Shalom, & Dayan, 1998). Determining the solutions of the algebraic equation and the Riccati difference equation when the system can be observed or detected but not necessarily stabilized in the filter is discussed by Siew Wah Chan et al. in (Chan et al., 2002), C.E DeSouza et al. have done a further analysis in (De Souza, Gevers, & Goodwin, 1986) for the case of singular state transition matrices. Upper and lower bounds can be used to evaluate the performance of a tracking system or as a means for the scheduling of sensors in real time. Derivation of the error covariance upper and lower bounds along with the criterion for the exponential by varying IMM algorithm's stability for one of the classes of the Markov jump linear system is already presented (Hwang, Seah, & Lee, 2017). In the same manner, the asymptotic stability and exponential stability of Riccati equation are analyzed in (Agarwal, 1992;Deyst & Price, 1968) and (Wang & Guo, 1999) respectively. Zhou Fucheng (Fucheng, 2010) explains the undecimated discrete wavelet transformation of the lifting system to analyze wind turbine gearbox failure diagnosis. A strong analysis of IMMKF for target tracking with constant acceleration and constant jerk model is performed by V.P.S Naidu (Naidu, 2009(Naidu, , 2010Naidu & Raol, 2008). Diagnosis of the fault signal is done by Zhou Wen-jing (Wen-jing, Yan-xia, & Long, 2012) and Zijun Zhang (Zhang, Verma, & Kusiak, 2012) with the help of wavelet analysis of time-frequency characteristics. A unique case of the IMM algorithm with a simple mixing probability, is the IMM algorithm to which the stability analysis is reasonably applied, is presented in this article. The behavior of the upper and lower bound performance of the IMM Kalman Filter is demonstrated by virtue of two examples namely (i) a maneuvering target tracking problem in air traffic control with constant acceleration,constant jerk models and (ii) a defective wind turbine gear box. One can refer (Bar-Shalom & Rong Li, 1995) to understand the basic principles and techniques on tracking. For further study on the literature the reader can refer (Anderson & Moore, 2005;Bar-Shalom, Rong Li, & Kirubarajan, 2001;Li & Bar-Shalom, 1993;Mehrotra & Mahapatra, 1997;Sinopoli et al., 2004). To the best of our knowledge, the investigation of jerk data obtained from the vibration of acceleration data via the IMM Kalman filter is a novel endeavor in the model monitoring to diagnose changes in the vibration excitation of the gearbox.

Motivation and contribution
In previous literature, most of the researchers use AE 2 ffiffiffiffiffiffiffiffiffiffi f P -ð k j kÞ q as a theoretical bounds to verify state error in (Naidu, 2009). In other words, the bounds of state error ð˘À^Þ is computed with the help of state error covariance matrix. When the state errors are within the bounds, then the uncertainty in the state estimates is limited and the filter is robust. This paper focuses on analyzing the stability of IMMKF algorithm by verifying the upper and lower bounds of state error covariance in (Hwang et al., 2017).The stability of IMMKF algorithm through upper and lower bounds of state error covariance matrix is determined using the Equations (21)-(33) .It is observed that values of the state error covariance are within the bounds, which shows that the IMMKF is more stable and robust. In this paper, at the beginning, the major contribution of the above-mentioned algorithm using IMMKF has been applied for tracking maneuvering target then it is extended by tracking vibration of wind turbine gearbox. At this moment, to meet the increasing energy demand in the electric power sector, renewable energy system plays a leading role in power system. One of the major sources of renewable energy is wind energy. In wind energy system, gearbox acts as a main component. The fault signals of the gearbox are non-stationary. The vibration of wind turbine gearbox has to be tracked accurately in order to diagnose the fault in gearbox of wind turbine.
The current research article is framed as follows. In section 2, we introduce the stability of IMM algorithm through upper and lower bound while the examination of a maneuvering target using the Constant Acceleration Model and Constant Jerk Models is presented in section 3. Identification of the failed components by the usage of the IMMKF algorithm with the upper and lower bounds of the state error covariance is investigated in section 4. The computational requirements of the algorithm are presented in section 5. Conclusions are elucidated in Section 6. It is observed that the suggested heterogeneous step algorithms with upper and lower bounds are more desirable for tracking maneuvering target and for failed component identification.

The IMMKF algorithm
The IMM algorithm incorporates r interacting filters (see Figure 8) functioning simultaneously with every filter corresponding to a model of the following stochastic hybrid system, (Hwang et al., 2017) where˘2 < n is the state of every system and z 2 < P is the measurement vector.
F mðkÞ and H mðkÞ are the matrices of the system corresponding to every model mðkÞ 2 f1; 2 . . . rg at time k : w mðkÞ ð kÞ and v mðkÞ ð kÞ are uncorrelated Gaussian noise vectors with white zero mean, Q mðkÞ and R mðkÞ are the respective covariance matrices. The evolvement of the model mðkÞ is denoted by ½ symbolize a conditional probability. In the above system, for all _ a; where I is the identity matrix.
Let Z k :¼ zð1Þ; zð2Þ; Á Á Á zðkÞ f g be the set of measurements up to time k. The IMM algorithm calculates the relative posterior mean _ _ b ð kÞ and P 2 _ b ðkÞ for each Kalman filter _ b, and the mode recursively as follows in (Agarwal, 1992). Let us suppose that, . . r are estimated from the last iteration at time k À 1 then the IMM algorithm computes the following for every k: • Interaction/Mixing Calculate the mixing probability For every Kalman filter _ b the initial conditions are calculated as denotes the residual covariance, and P 2 _ b ð kÞ; ðP 2 _ b ð kj k À 1ÞÞ is the posterior (prior) state covariance.

• Model Probability Update
Herewith Þis a p-dimensional multivariate Gaussian probability density function along with mean zero and ∑ covariance. The model probability is then given by • State Estimate and Covariance combiner
Proof. By the replacement of (18) and (19) into (9), we have with a change in notations of the The function Υ From (21) and (22), where By analyzing (16) and (24), we see that P 2 Ã ð kÞ can be considered as the Kalman filter covariance of the system(12)-(13) with Fð kÞ ¼ F _ b k ; Fð k À 1Þ ¼ F _ b kÀ1 , and so on. Hence by the above Lemma 1, where β _ a is the positive scalar in a such way that where _ a is the dependency on a given model sequence (the _ a-th desirable realization of the Thus, P 2 _ b k ð kÞ ! β 1 I where β 1 ; minfβ 1 ; β 2 ; Á Á Á β r k g.

An upper bound for
From (7) the upper bound of the covariance P 2 _ b ð kÞ is derived where K _ b ð kÞ of the IMM algorithm is given by Furthermore, We get We define a function χ _ a; k ðP 2 ; Theorem 2. If the system (1)-(2) is uniformly controllable and uniformly observable, then,for all By iterations, we have (with a change in the subscript notations) T Using Eð kÞ Eð k À 1Þ Á Á Á Eð k 0 Þ c 0 ς kÀ k 0 0 , c 0 > 0; 0 0 1 and considering By taking 2-norm on (33), and utilizing ∑ r (34), we see  (29) and (9) we see that P 2 _ b ð kÞ is uniformly bounded from above.  Due to these transformations, the elements of measurement noise covariance in the Cartesian coordinates become correlated. The elements in noise covariance matrix can be updated using the following relations where ς 2 r ; ς 2 # and ς 2 Ψ denote the variances of the measurements noise in r; # and Ψ dimensions, respectively.

Acceleration model
Let us define the acceleration model as X ac ð k þ 1Þ ¼ F ac X ac ð kÞ þ ω ac ð kÞ; Z ac ð k þ 1Þ ¼ H ac X ac ð kÞ þ υ ac ð k þ 1Þ: Also H ac ¼ The state transition and noise covariance matrices for the acceleration model are

Derivation of H ac
Calculation of the Jacobian measurement matrix H ac is obtained by where i ¼ 1; 2 Á Á Á is the length of the measurement vector Z ac and j ¼ 1; 2 Á Á Á length of the state vector X ac

Jerk model
Let us define the jerk model as where u 1je ; u 2je ; Á Á Á are the driving measurement noise for the jerk model. For the jerk model, we illustrate F je as the state transition matrix and Q je as noise covariance matrix such as where

Wind turbine
One of the major components in any wind turbines is the gearbox as shown in the Figure 1. The fault signals of the gearbox are commonly non-stationary and extremely polluted with noise. To validate the gearbox fault diagnosis approach, vibration signals collected from a real wind turbine gearbox are analyzed. The experimental vibration data are collected from the defective gearbox of wind turbine under test. This enables in the determination of the key factor for the defective gear box and to develop the monitoring strategies. Accelerometers mounted on 12 locations exterior of the gearbox, generator, and main bearing collected the vibration acceleration with 40 kHz sampling rate. The 3-D coordinate system describes the direction of the drivetrain vibration acceleration and is observed by accelerometers. Figure 2 illustrates the vibration acceleration coordinate system. The point of intersection of plane and the planet-carrier rotation dividing the torque arm cylinder midway along their length is the origin of the coordinate system. The x axis describes the acceleration of the system along the main shaft axis and the downwind side, with the y-axis reflecting the direction of vibration acceleration, horizontally perpendicular to the x-axis. The x,y and z-axis are orthogonal to each other. Although the vibration acceleration of the system is depicted by a 3-D coordinate system, the mounted accelerometers can only sense a maximum of two acceleration directions. Investigation of the failed components and determination of the root causes of failure are performed by the time and frequency-domain analysis of the vibration acceleration data. Conversion of the vibration, the high-frequency jerk data (40kHz), into much lower frequency data 1 15 Hz À Á is achieved by the computation of the value of mean jerk at 15s intervals along with the standard deviation and the maximum value of the jerk data. The chosen sampling frequency of 1 15 Hz, reduces the data size for ease in time-domain analysis.

Case 1
The maneuvering target motion trajectory has been simulated using CA, CJ models for a period of 250s with sampling time T as 1sec. The motion of the target begins at 0; 0; 100 ð Þwith a constant velocity of 100 m/s and 10 m/s along the x-and y-axes. A step jerk of 0.009m/s 3 is applied along the x-axis at 50s and −0.009m/s 3 is applied along z-axis at 125s, resulting in a 1. cubic position 2 parabola velocity, and 3.ramp acceleration variation. The following variances are added to the range, azimuth, and elevation of the random noises to generate noise measurement data.
Both models have a correlation factor α that is chosen as 0.006 to maintain a degree of analytical similarity.Correlation of the target acceleration model is a requirement for acceleration model, while the jerk model has the target jerk correlated. The acceleration model has its process noise variance Q ac ¼ 2ασ 2 ac ; with σ ac ¼ 1:8m=s 2 where as for the jerk model, process noise variance Q je ¼ 2ασ 2 je , with σ je ¼ 0:009m=s 3 are used. The simulated trajectory of position,velocity,acceleration and jerk along x,y,z axes is shown in Figure 3. The average of 75 Monte Carlo simulations was utilized to analyze the performances of both the models. Figure 4 depicts the estimated and true trajectories. Figure 5

Case 2
The jerk data are measured by the sensors AN9, which is located in high speed shaft downwind bearing radial with acceleration direction in the positive direction of the Z axis with the period of 10-min test. The initial state starts at jerk (0, 0, 3.55) which is measured by AN9 from Figure 3 in (Li & Jilkov, 2010). All through the test experiment, there is a progressive increase in the change of vibration. Generation of noisy measurement data is done by the addition of random noises with the following variances to the range, azimuth, and elevation given by Range ðς r Þ 2 : 2:25e À 1m 2 ; Azimuth ðς # Þ 2 : 2e À 4 À rad 2 ;  Elevation ðς Ψ Þ 2 : 2e À 4 À rad 2 .
As in the prior case, correlation factor α is chosen as 0.006 for the models. Process noise variance of Q ac ¼ 2ασ 2 ac ; σ ac ¼ 1:5m=s 2 is used for acceleration model while the jerk model, a process noise variance of Q je ¼ 2ασ 2 je ; σ je ¼ 0:0015m=s 3 is used. Figure 6 depicts the estimated and true trajectories. Figure 7

Summary and conclusions
A two-model IMM Kalman filter algorithm with constant acceleration and constant jerk is carried out in tracking maneuvering target in aircraft. The boundedness of Riccati equation is analyzed by determining the upper and lower bounds of state error covariance, which helps to decide the stability of the algorithm. The performance of the IMMKF algorithm is also evaluated for tracking the vibration obtained from the wind turbine gearbox. This helps us to diagnose the fault in the gearbox of wind turbine. Finding the upper bound and lower bound of the state error covariance is more effective in identifying the fault diagnosis. The proposed work can be extended by investigating various applications using extended or unscented Kalman filter. In addition, the algorithm can be modified for tracking of multiple objects through the properties of Riccati difference equation.