© AFM, EDP Sciences 2018

1 Introduction

Order tracking using the state space approach is one of the tools widespread for the processing of non-stationary signals. The so-called the state space model is composed of two equation: the state equation and the measurement equation. The technique the most presented in this area is the Kalman estimator and more precisely the Vold-kalman estimator in the area of the mechanical systems diagnosis [1]. Vold et al. present the theoretical basis about this estimator in [2]. This kind of estimator suppose that the measurement noise and the process noise are centered, Gaussian and white with known statistics.

In the literature many works on the Vold-kalman estimator for order tracking has been done. Pan and Lin have realized an interesting explorative study on the Vold-kalman estimator [3]. Behrouz and al. also applied this estimator to diagnose a bearing default and they have translated the state equation in term of second order autoregressive model [4]. These study have provided conclusive results. However, the unrealistic assumptions on the noises naturally limit the application of this estimator in real cases.

Therefore, the H_∞ estimator is proposed in this paper to evaluate the amplitude modulation and the phase modulation. To estimate these modulations the vibration signal is described using the state variables. Then these latter are modelled by a Taylor series. This method generalize that of Vold-kalman. With the H_∞ estimator we make no assumption on the noise statistics. They must only be of finite energy. More details on the discrete H_∞ estimator can be found in the work of Shen and Deng [5].

This paper is structured as follows: Section 2 presents the theoretical foundation about the H_∞ estimator and Section 3 provides an example of simulation which validated our proposal.

2 Theoretical background

2.1 Problem formulation

In this paper, the gearbox vibration signal is modelled as $$y\left(t\right)={\displaystyle {\displaystyle \sum }_{i=1}^M}{A}_i\left(t\right)\cos \left(2\pi {\displaystyle \underset{0}{\overset{t}{{\displaystyle \int }}}}{f}_i\left(u\right)du+{\phi }_i\left(t\right)\right)+v(t)$$(1) where A_i and φ_i are respectively the amplitude and the phase of the ith order, v is the measurement noise which contains the unwanted part of the signal, f_i = o_if_r is the instantaneous frequency of the order of interest with f_r the reference frequency and o_i the value of the orderi.

In the discrete form, (Eq. (1)) becomes: $$y\left(k\right)={\displaystyle \sum }_{i=1}^M{A}_i\left(k\right)\cos \left({\theta }_i\left(k\right)+{\varphi }_i\left(k\right)\right)+v(k),k=0,2,\dots ,n-1,$$(2) where ${\theta }_i\left(k\right)=2\pi {\displaystyle {\displaystyle \sum }_{j=1}^k}\frac{f_i(j)}{f_s}$ is the angular displacement and f_s is the sampling frequency.

The purpose is to estimate the amplitude and the phase of some specific orders of interest using the H_∞ estimation approach. For this, the problem is formulated in term of estimation of the state variables. Note that the amplitude and the phase modulation are the key features to diagnose the gear state [6].

2.2 State space modelling

Let us consider the formula established in (Eq. (2)). The purpose here is to build the measurement and the state equation.

Linearizing (Eq. (2)) leads to: $$y\left(k\right)={\displaystyle {\displaystyle \sum }_{i=1}^M}\left[\begin{array}{cc}\hfill \cos ({\theta }_i\left(k\right))\hfill & \hfill -\sin ({\theta }_i\left(k\right))\hfill \end{array}\right]\left[\begin{array}{c}\hfill a_{i,c}(k)\hfill \\ \hfill a_{i,s}(k)\hfill \end{array}\right]+v(k)\text{,}$$(3)

where a_i,c = A_icosφ_i and a_i,s = A_isinφ_i. Let put $a_i(k)=\left[\begin{array}{c}\hfill a_{i,c}(k)\hfill \\ \hfill a_{i,s}(k)\hfill \end{array}\right]$ and $B_i\left(k\right)=\left[\begin{array}{cc}\hfill \cos \left({\theta }_i(k)\right)\hfill & \hfill -\sin \left({\theta }_i(k)\right)\hfill \end{array}\right]\text{.}$

The amplitudes a_i,c and a_i,s are unknown. For estimating them, these amplitudes are modeled by a polynomial approximation as follows: $$a_{i,c}(k)={\displaystyle {\displaystyle \sum }_{q=0}^N}{\alpha }_{i,c}^q(k)t^q(k)\text{,}$$(4) $$a_{i,s}(k)={\displaystyle {\displaystyle \sum }_{q=0}^N}{\alpha }_{i,s}^q(k)t^q(k),\mathrm{i}=1,2,\dots ,\mathit{M}$$(5) and the coefficients of the polynomial by a random walk process such as: $${\alpha }_{i,c}^q(k+1)={\alpha }_{i,c}^q(k)+w_{i,c}(k)\text{,}$$(6) $${\alpha }_{i,s}^q(k+1)={\alpha }_{i,s}^q(k)+w_{i,s}(k)\text{,}$$(7)where w_i,. is a random signal. With those new variables (3) can be rewritten as: $$y(k)={\displaystyle {\displaystyle \sum }_{i=1}^M}{B}_i(k){\displaystyle \tilde{T}}(k)x_i(k)+v(k)\text{,}$$(8)with ${\displaystyle \tilde{T}}(k)={\left[\begin{array}{cc}\hfill T(k)\hfill & \hfill T(k)\hfill \end{array}\right]}^T\text{,}$ $T(k)=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill t(k)\hfill & \hfill \cdots \hfill & \hfill t^N(k)\hfill \end{array}\right]$ and $x_i\left(k\right)=\left[\begin{array}{c}\hfill x_{i,c}\hfill \\ \hfill x_{i,s}\hfill \end{array}\right]$ with $x_{i,\mathrm{.}}={\left[\begin{array}{cccc}\hfill {\alpha }_{i,\mathrm{.}}^0\hfill & \hfill {\alpha }_{i,\mathrm{.}}^1\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{i,\mathrm{.}}^N\hfill \end{array}\right]}^T\text{.}$Note that A^T is the transpose of the matrix A.

Assuming that the measurement matrix is $H(k)=\left[\begin{array}{cccc}\hfill B_1(k){\displaystyle \stackrel{ }{T}}(k)\hfill & \hfill B_2(k){\displaystyle \stackrel{ }{T}}(k)\hfill & \hfill \cdots \hfill & \hfill B_M(k){\displaystyle \stackrel{ }{T}}(k)\hfill \end{array}\right]\text{,}$ the following measurement equation is obtained: $$y(k)=H(k)x(k)+v(k)\text{,}$$(9) where $x(k)={\left[\begin{array}{cccc}\hfill x_1(k)\hfill & \hfill x_2(k)\hfill & \hfill \cdots \hfill & \hfill x_M(k)\hfill \end{array}\right]}^T$ and v is the measurement noise with a covariance matrix V.

Then the state equation is: $$x(k+1)=Fx(k)+w(k)\text{,}$$(10) where $F=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill \end{array}\right]$ and w(k) = [w_1,c w_1,s w_2,c w_2,s ⋯ ⋯ w_M,c w_M,s]^T is the process noise with a covariance matrix W.

2.3 Discrete H_∞ estimator design

From equations (9) and (10) the following state space model is considered: $$\{\begin{array}{c}\hfill \begin{array}{cc}\hfill x_{k+1}=\hfill & \hfill F{x}_k+B{w}_k\hfill \end{array}\hfill \\ \hfill \begin{array}{cc}\hfill y_{k+1}=\hfill & \hfill H_k{x}_k+v_k\hfill \end{array}\hfill \end{array}\text{.}$$(11) Let us note $e_k=x_k-{{\displaystyle \hat{x}}}_k$ the estimation error where ${{\displaystyle \hat{x}}}_k$ is the estimate of x_k and E{ . } will stand for the expectation value.

Several facts may be used against the Kalman estimator although it is an attractive and powerful tool to estimate x_k;

• the Kalman estimator minimizes $E\left\{e_k{e}_k^T\right\}$ while the user may be interested in minimizing the worst-case error;

• the Kalman estimator assumes that the noises are zero-mean with Gaussian distribution;

• the Kalman estimator assumes also that $E\left\{v_k{v}_k^T\right\}$ and $E\left\{w_k{w}_k^T\right\}$ are known.

These limitations have led to the statement of the H_∞ estimation problem. Several formulations exist in the literature. The H_∞ estimator solution that we present here is originally developed by Ravi Banar [7] and further explored by Shen and Deng [5]. These pioneers define the following cost function: $$J=\frac{{{\displaystyle \sum }}_{k=0}^{n-1}{\Vert x_k-{{\displaystyle \hat{x}}}_k\Vert }_{\mathit{Q}}^2}{{\Vert x_0-{{\displaystyle \hat{x}}}_0\Vert }_{{\mathit{p}}_0^{-1}}^2+{{\displaystyle \sum }}_{k=0}^{n-1}\left({\Vert {\mathit{w}}_{\mathit{k}}\Vert }_{{\mathit{W}}^{-1}}^2+{\Vert {\mathit{v}}_{\mathit{k}}\Vert }_{{\mathit{V}}^{-1}}^2\right)}\text{,}$$(12) where ${{\displaystyle \hat{x}}}_0$ is an estimate of x₀, Q > 0, P₀ > 0, W > 0 and V > 0 are the weighting matrices and are left to the choice of the designers and depend on the performance requirements. The notation ${\Vert {\mathit{x}}_{\mathit{k}}\Vert }_{\mathit{Q}}^2$ defines the weighted Q − L₂ norm, i.e, ${\Vert {\mathit{x}}_{\mathit{k}}\Vert }_{\mathit{Q}}^2=x_k^{T}Q{x}_k\text{.}$

Problem statement [8]: Given the scalar γ > 0, find estimation strategy that achieve $$\sup J<1/\gamma$$(13) where “ sup” is the supremum value and γ is the desired level of noise attenuation.“

The H_∞ estimation problem consists of the minimization of the worst possible amplification of the estimation error. This can be interpreted as a “minmax” problem in which the estimation error is to be minimized and the exogenous disturbances (w_k and v_k) and the error of initialization ($x_0-{{\displaystyle \hat{x}}}_0$) is to be maximized.

Remember that unlike the Wiener/Kalman estimator, the H_∞ estimator deals with deterministic noises and no a priori information on their statistic properties are required. The solution of the H_∞ estimation problem is given in the theorem below from [5].

Theorem: Let γ > 0 be a prescribed level of noise attenuation. Then, there exists a H_∞ estimator for x_k if and only if there exists a stabilizing symmetric solution P_k > 0 to the following discrete-time Riccati equation: $$P_{k+1}=F{P}_k{\left(I-\gamma Q{P}_k+H_k^T{V}^{-1}{H}_k{P}_k\right)}^{-1}{F}^T+BW{B}^T\text{.}$$(14) Then the H_∞ estimator gives the estimate ${{\displaystyle \hat{x}}}_k$ of x_k such as: $${{\displaystyle \hat{x}}}_{k+1}=F{{\displaystyle \hat{x}}}_k+K_k\left(y_k-H_k{{\displaystyle \hat{x}}}_k\right),{{\displaystyle \hat{x}}}_0=x_0\text{.}$$(15)K_k is the gain of the H_∞ estimator and is given by: $$K_k=F{P}_k{\left(I-\gamma Q{P}_k+H_k^T{V}^{-1}{H}_k{P}_k\right)}^{-1}{H}_k^T{V}^{-1}\text{.}$$(16)

Another way to solve the Riccati equation (14) is presented by Yaesh and Shaked [9]. The method is given as follows:

• 1.

  Form the Hamiltonian

  $$Z=\left[\begin{array}{cc}\hfill F^{-T}\hfill & \hfill F^{-T}[H^T{R}^{-1}H-\gamma I]\hfill \\ \hfill BQ{B}^T{F}^{-T}\hfill & \hfill F+BQ{B}^T{F}^{-T}[H^T{R}^{-1}H-\gamma I]\hfill \end{array}\right]\in {ℝ}^{2n*2n}\text{,}$$(17)

  where  n  is  x  dimension.

• 2.

    Find the eigenvectors of Z corresponding to the eigenvalues  ℰ_i (i = 1, … , n)    outside the unit circle

• 3.

    Form the matrix of the corresponding eigenvectors denoted by:

  $$({\epsilon }_1{\epsilon }_2\cdots {\epsilon }_n)\equiv \left[\begin{array}{c}\hfill {\mathcal{X}}_1\hfill \\ \hfill {\mathcal{X}}_2\hfill \end{array}\right],{\mathcal{X}}_1{\mathcal{X}}_2\in {\mathrm{ℝ}}^{n*n}\mathrm{.}$$(18)

• 4.

  Compute $P={\mathcal{X}}_2{\mathcal{X}}_1^{-1}\text{.}$

Note that the smaller γ, the more easy the problem is to solve. When γ tends to ${\mathit{\gamma }}_{\mathrm{o}\mathrm{p}\mathrm{t}}$ (the optimal value of  γ) the eigenvalues of P tend to infinity and therefore ${\mathcal{X}}_1$ is close to a singular matrix. Shaked and Theodor [10] investigated the behavior of the optimal H_∞ estimator when γ tends to ${\mathit{\gamma }}_{\mathrm{o}\mathrm{p}\mathrm{t}}$. They showed that when  γ  reaches ${\mathit{\gamma }}_{\mathrm{o}\mathrm{p}\mathrm{t}}$, there exists at least one or more unbounded eigenvalues.

In the special case, where γ → 0, the H_∞ estimator reduces to a Kalman estimator.

3 Numerical implementation

In this section, a synthetic signal is used to illustrate the performances of H_∞ estimation approach. The generated signal (see Fig. 1) is described by the following equation. $$y(t)={\displaystyle {\displaystyle \sum }_{i=1}^3}{A}_i(t)\text{cos}\left(2\pi o_i{\displaystyle \underset{0}{\overset{t}{{\displaystyle \int }}}}{f}_r(u)du\right)+v(t)\text{,}$$(19) where f_r is the instantaneous frequency linearly increasing from 0 to 50  Hz in 5 s, o_i contains the order's number and v is the measurement noise. The signal is composed of three orders presented in the Table 1. Figure 2 displays the rpm-frequency spectrum using the conventional windowing Fourier transform that characterizes three orders.

The results presented below have been got using a Monte-Carlo simulation based on 400 iterations.

The parameters of the estimator have been taken as follows:

• the covariance of the process noise W = 10⁻⁹;

• the covariance of the measurement noise V = 10⁻³;

• the initial covariance error P₀ = 10⁻³;

•  the level of the noise attenuation γ = γ_opt = 10^0.178.

γ_opt is equal to the greatest value that guarantees the stability of the matrix P. This stability is reached, according to Yaesh and Shaked [11], when Ps eigenvalues are bounded in the unit circle. As plotted in the Figure 3, this stability is reached for γ = 10^0.178. Beyond this value there exists at least one or more eigenvalues that are outside the unit circle.

The measurement noise is modelled by a Poisson noise as mentioned in [11]. The Kalman estimator algorithm presented by Dan Simon [12] and the H_∞ estimator have been applied to the generated vibration signal. The performance of both estimators is measured in term of signal to noise ratio. Table 2 gives the performance got for the two estimators. In both cases the H_∞ estimator provides a better result than the Kalman estimator. The SNR_out value is the signal to noise ratio calculated by: $$SN{R}_{\text{out}}=10*{\log }_{10}\frac{{{\displaystyle \sum }}_{k=1}^N{y}_k^2}{{{\displaystyle \sum }}_{k=1}^N(y_k-{{\displaystyle \hat{y}}}_k{)}^2}\text{,}$$(20) where N is the number of samples, y_k is the noiseless signal at times k and ${{\displaystyle \hat{y}}}_k$ is the estimated or filtered signal. The criterion of comparison is improved by about 0.7  dB using the H_∞ estimator. Therefore the H_∞ estimator is a good alternative to deal with real situation where the noises are not really Gaussian.  

Figures 4–6 show the effectiveness of the H_∞ estimator for order tracking in non-stationary signal processing. We see in this last figure that the estimated we got by the H_∞ estimation is closer to the original amplitude than the Kalman estimation.

Fig. 1   Synthetic signal.

Fig. 2 Illustration of rpm-frequency spectrum.

Fig. 3 Maximum of the eigenvalues of the covariance matrix error.

Fig. 4 Amplitude of the 1st order estimated using the  H∞  and the Kalman estimator.

Fig. 5 Amplitude of 3rd order estimated using the  H∞  and the Kalman estimator.

Fig. 6 Amplitude of the 9th order estimated using the  H∞  and the Kalman estimator.

4 Conclusion

Through this paper a method has been developed to estimate order's amplitude based on the H_∞ estimation in non-stationary operations. This method uses the information of the instantaneous frequency of the signal and makes no assumption on the noises statistics. It takes advantage on the classical Kalman estimation and it can be consider as an extension of this last one. Since the estimator is designed to minimize the worst case-disturbances, the H_∞ estimation approach is more robust to process any kind of noisy signal. The application of this method in real-life data will concern our future research.

References

All Tables

All Figures

	Fig. 1   Synthetic signal.
In the text

	Fig. 2 Illustration of rpm-frequency spectrum.
In the text

	Fig. 3 Maximum of the eigenvalues of the covariance matrix error.
In the text

	Fig. 4 Amplitude of the 1st order estimated using the  H∞  and the Kalman estimator.
In the text

	Fig. 5 Amplitude of 3rd order estimated using the  H∞  and the Kalman estimator.
In the text

	Fig. 6 Amplitude of the 9th order estimated using the  H∞  and the Kalman estimator.
In the text