Failure prognosis of multiple uncertainty system based on Kalman filter and its application to aircraft fuel system

Abstract

Failure prognosis is the key point of prognostic and health management or condition-based maintenance, the multiple uncertainty sources in real world will lead to inaccurate prediction. In this paper, an advanced failure prognosis method with Kalman filter is presented to address the real-world uncertainties. The multiple uncertainty sources are analyzed and classified first and then theoretical methods are derived, respectively, for the different uncertainty sources. Afterward, the failure prognosis algorithm is developed by taking into consideration. In the end, an aircraft fuel feeding system health monitoring case simulation is presented to demonstrate the effectiveness of the proposed method.

Introduction

Various failure prognosis methods based on measured data or degradation information can be used to evaluate system’s degradation trend and estimate the remaining useful life. Failure prognosis methods can be seen as the bridge between the front-end field sensor data acquisition and back-end fusion/inference machine and also the critical component of prognostic and health management (PHM) or condition-based maintenance (CBM).^1,2 High credible failure prognosis plays an important role in maintenance of interval planning, operation schedule, logistical support optimization, and so on. Moreover, for the safety-critical system such as aircraft or electrical system, it is essential to estimate the failure state exactly to avoid catastrophic accidents and even human lives losing.

Over the years, research on failure prognosis in different fields has been a topic of considerable interest and obtained great achievements.^3,4 Traditionally, the corresponding approaches have been classified into two categories, for example, model based and data driven. The former utilizes physical knowledge or failure mechanism of system to realize prognosis, which is based on system’s physics modeling or expensive accelerated life test, while the latter only relies on available measured data and statistical models.

However, no matter what kinds of method, if a large number of errors or uncertainties are contained during the prognosis process, the credibility of failure prognosis information will be suspected. To address these problems, it needs to consider the different uncertainty sources during every step in PHM’s realization and void the introduction of any “noise” that leads to inaccurate prediction.^2,5,6 Gu et al. ⁷ studied various uncertainty sources in the whole process and classified them into four kinds, for example, parameter, measurement, failure criteria, and future usage uncertainty, and with sensitivity analysis, they found that measurement items are the main source. Apart from these, Sun et al. ² grouped these uncertainties into three different categories: model uncertainty, measurement uncertainty, and uncertainties of the product’s characteristic parameters. For these problems, Gu et al. ⁸ developed a method for monitoring, recording, and analyzing the life cycle of vibration loads for remaining life prediction of electronics. Wang et al. ⁹ presented a degradation model for the prediction of the residual life based on the adapted Brownian motion-based approach with a drifting parameter.

The stochastic filtering–based method utilizes the current data’s new information to correct the estimation results.^10–12 If the measured data have the trend to be used to predict, then this method can be fit for failure prognosis in uncertain or unknown operating conditions. The Kalman filter is a well-known and widely used stochastic filtering–based method which was proposed by R.E. Kalman ¹³ as early as the 1960s; it is a linear unbiased and minimum variance estimator under the linear assumption. However, most practical systems are nonlinear, and later, a variety of methods were developed to deal with nonlinear filtering problems, such as extended Kalman filter (EKF), unscented Kalman filter (UKF), central difference Kalman filter, and particle filter, ¹⁴ in which EKF is a generalized and main form of Kalman filter for nonlinear system, and it realize filtering calculation through system linearization and higher-order item omission, which can solve the deviation or uncertainty problem between actual system and nominal system. And the other way to solve such problem is considering the multiple uncertainties for failure prognosis engineering practice, which is also the origin of this article in theoretical aspect, it is necessary to consider all the uncertainties in the whole process and introduce some methods to reduce their impacts on prognosis.

The outline of this article is as follows: the failure prognosis techniques based on filtering algorithm are presented first in section “Fault prognosis theory based on stochastic filtering” and then in section “Uncertainty source analysis,” the different uncertainty sources and their impacts on fault prognosis are analyzed. In section “Failure prognoses with multiple uncertainty sources,” several algorithms are conducted to reduce the impact of uncertainties. Furthermore, an integrated failure prognosis algorithm is carried out in section “Failure prognosis value estimation algorithms.” A case simulation about aircraft fuel system is provided in section “Simulation examples” to show the validity of the proposed approach. In the end, section “Conclusion” presents several comments and final remarks.

Fault prognosis theory based on stochastic filtering

Regardless of which methods we choose for failure prognosis, the measured data must have the predictable trend. The stochastic filtering–based method updates the state estimation based on the new information from newest measured sensor data and then the future state value is predicted by statistical tracking algorithms with the assumption of no measurement noise. In this section, the classical failure prognosis methods based on standard Kalman filter is presented first ¹⁵ without uncertainty.

Kalman filtering

Consider the discrete state–space model equation as^15–17

$$\{\begin{array}{c}\hfill {\mathit{x}}_{\mathit{k}}={\mathit{A}}_{\mathit{k}-1}{\mathit{x}}_{\mathit{k}-1}+{\mathit{B}}_{\mathit{k}-1}{\mathit{v}}_{\mathit{k}-1},\mathit{k}\ge 0\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_{\mathit{k}}={\mathit{C}}_{\mathit{k}}{\mathit{x}}_{\mathit{k}}+{\mathit{w}}_{\mathit{k}}}\end{array}$$

(1)

where ${\mathit{A}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$ is the system matrix, ${\mathit{y}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{m}}$ is the measurement output, ${\mathit{B}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{n}\times \mathit{r}}$ is the noise matrix, ${\mathit{x}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{n}}$ is the state vector, ${\mathit{C}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{m}\times \mathit{n}}$ is the measurement matrix, ${\mathit{v}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{n}}$ is the system noise, and ${\mathit{w}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{m}}$ is the measurement noise.

The system noise ${\mathit{v}}_{\mathit{k}}$ and measurement noise ${\mathit{w}}_{\mathit{k}}$ satisfy ¹⁶

$$\begin{array}{c}\hfill \mathit{E}\left[{\mathit{v}}_{\mathit{k}}\right]=0,\mathit{Cov}[{\mathit{v}}_{\mathit{k}},{\mathit{v}}_j]=\mathit{E}\left[{\mathit{v}}_{\mathit{k}}{\mathit{v}}_j^{\mathit{T}}\right]={\mathit{Q}}_{\mathit{k}}\hfill \\ \multicolumn{1}{c}{\mathit{E}\left[{\mathit{w}}_{\mathit{k}}\right]=0,\mathit{Cov}[{\mathit{w}}_{\mathit{k}},{\mathit{w}}_j]=\mathit{E}\left[{\mathit{w}}_{\mathit{k}}{\mathit{w}}_j^{\mathit{T}}\right]={\mathit{R}}_{\mathit{k}}}\\ \multicolumn{1}{c}{\mathit{Cov}[{\mathit{v}}_{\mathit{k}},{\mathit{w}}_j]=\mathit{E}\left[{\mathit{v}}_{\mathit{k}}{\mathit{w}}_j^{\mathit{T}}\right]=0}\end{array}\}$$

(2)

where ${\mathit{Q}}_{\mathit{k}}$ is the system noise’s variance matrix and assumed to be a nonnegative matrix and ${\mathit{R}}_{\mathit{k}}$ is the measurement noise’s variance matrix and assumed to be a positive matrix.

For given measurement output ${\mathit{y}}_{\mathit{k}}(\mathit{k}=0,1,\dots ,\mathit{n})$ , the optimal estimation ${\hat{\mathit{x}}}_{\mathit{k}}$ for system state ${\mathit{x}}_{\mathit{k}}$ can be calculated to satisfy the object function ¹⁶

$$\mathit{J}=\mathit{E}\left\{{\tilde{\mathit{x}}}_{\mathit{k}}{\tilde{\mathit{x}}}_{\mathit{k}}^{\mathit{T}}\right\}=\mathit{E}\left\{[{\mathit{x}}_{\mathit{k}}-{\hat{\mathit{x}}}_{\mathit{k}}]{[{\mathit{x}}_{\mathit{k}}-{\hat{\mathit{x}}}_{\mathit{k}}]}^{\mathit{T}}\right\}$$

(3)

And then the standard Kalman filter equation can be obtained as^16,17

$$\begin{array}{l}{\widehat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}={\mathit{A}}_{\mathit{k}-1}{\widehat{\mathit{x}}}_{\mathit{k}-1}\\ {\widehat{\mathit{x}}}_{\mathit{k}}={\widehat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}\left({\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}{\widehat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}\right)\\ {\mathit{K}}_{\mathit{k}}={\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\mathit{C}}_{\mathit{k}}^{\mathit{T}}{\left({\mathit{C}}_{\mathit{k}}{\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\mathit{C}}_{\mathit{k}}^{\mathit{T}}+{\mathit{R}}_{\mathit{k}}\right)}^{-1}\\ {\mathit{P}}_{\mathit{k}/\mathit{k}-1}={\mathit{A}}_{\mathit{k}-1}{\mathit{P}}_{\mathit{k}-1}{\mathit{A}}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{B}}_{\mathit{k}-1}{\mathit{Q}}_{\mathit{k}-1}{\mathit{B}}_{\mathit{k}-1}^{\mathit{T}}\\ {\mathit{P}}_{\mathit{k}}=\left(I-{\mathit{K}}_{\mathit{k}}{\mathit{C}}_{\mathit{k}}\right){\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\left(\mathit{I}-{\mathit{K}}_{\mathit{k}}{\mathit{C}}_{\mathit{k}}\right)}^{\mathit{T}}+{\mathit{K}}_{\mathit{k}}{\mathit{R}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}^{\mathit{T}}\end{array}$$

(4)

where ${\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}$ is the state prediction in one step, ${\hat{\mathit{x}}}_{\mathit{k}}$ is the state estimation, ${\mathit{K}}_{\mathit{k}}$ is Kalman filtering gain, ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ is the mean-squared error of state prediction in one step, and ${\mathit{P}}_{\mathit{k}}$ is the mean-squared error of state estimation. With the predefined initial state ${\mathit{x}}_0$ and mean-squared error of state estimation ${\mathit{P}}_0$ , the filtering algorithm can be completed step by step as shown in Figure 1, which shows that the Kalman filtering method can predict system state one time step ahead based on the new information from the newest measurement ${\mathit{y}}_{\mathit{k}}$ .

Figure 1.

Flowchart of Kalman filtering.

Failure prognosis approach

Theoretically, all the failure prognoses are probability estimation method, which means that compared with current time, the further the failure prognosis, the greater the variance, and the worse the accuracy. Figure 2 illustrates this viewpoint well, where the x-axis is the system operation time (note: only for illustration) and the y-axis is the probability density function (PDF) of system failure during the whole operation time. However, as the possible fault point is close to the current time, failure prognosis error will decrease, and the failure prognosis value will converge to the actual failure time value, and PDF curve will become steeper and steeper, which means that the possibility of failure is getting higher and higher.

Figure 2.

Failure prognosis value’s probability density function in different time lengths.

The above Kalman filter utilizes the system motion model and new information in current time to estimate system’s state in the next time step. So, we can use the current estimated states, assume that the noise is zero and that there is no measured new information, and then extrapolate or predict step by step. Take one-step prediction as an example

$$\begin{array}{c}\mathrm{State\; estimation\; or\; extrapolation}:{\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}={\mathit{A}}_{\mathit{k}-1}{\hat{\mathit{x}}}_{\mathit{k}-1}\hfill \\ \mathrm{Covariance\; estimation}:\hfill \\ {\mathit{P}}_{\mathit{k}/\mathit{k}-1}={\mathit{A}}_{\mathit{k}-1}{\mathit{P}}_{\mathit{k}-1}{\mathit{A}}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{B}}_{\mathit{k}-1}{\mathit{Q}}_{\mathit{k}-1}{\mathit{B}}_{\mathit{k}-1}^{\mathit{T}}\hfill \end{array}$$

(5)

Such calculation will be iterated until the estimated value is equal to or beyond the predefined characteristic parameter threshold. Moreover, the mean and variance of estimation value are defined as the state and covariance estimation mentioned above, that is, at k th filtering step, the PDF of the estimated value is assumed to be

$$\mathit{f}\left({\mathit{x}}_{\mathit{k}}\right)=\frac{1}{\sqrt{2\pi {\mathit{p}}_{11}}}{\mathit{e}}^{\frac{-{({\mathit{x}}_{\mathit{k}}-{\hat{\mathit{x}}}_{\mathit{k}})}^2}{2{\mathit{P}}_{11}}}$$

(6)

where ${\mathit{p}}_{11}$ is the (1, 1) element of covariance matrix $\mathit{P}$ . With equation (6), the probability that system failed before this time point is shown as the shadow zone in Figure 3. Different from Figure 2, the x-axis of Figure 3 represents the estimated system state, and the y-axis indicates the failure PDF of the corresponding x-axis state value. In Figure 3, the first thick solid line in the left represents the current filtering value of the system, for example, in the right of this line, all the computing results are predicted value of the system state, and the second thick solid line in the right represents the k th step predicted value, and the corresponding solid curve represents PDF of the system failure in that time. Then, the failure probability can be calculated by the integral from the threshold point to positive infinity. Using this method, the system failure probability of different extrapolation points or estimated states before this time point can be calculated.

Figure 3.

Probability of the dedicated failure prognosis.

Uncertainty source analysis

According to the filtering estimation principle and flowchart in section “Kalman filtering,” the whole prediction process of stochastic filtering–based method can be summarized as shown in Figure 4. ¹⁸ The prediction process consists of the following steps:

1. Data sampling and transferring. Different types of sensors are used to collect corresponding system parameter. Meanwhile, the conversion of various physical signals to electrical signals is required in this step.

2. Data preprocessing and feature extraction. Data transmitted from sensors will be preprocessed to specified format for subsequent failure prognosis. Simplified, integrated, or compressed data will be the output in this step.

3. Determining the filtering model. As the different kinds of data have different characteristics, or the same data type in different system stages also has different characteristics, the filtering model of subsequent filtering calculation should be determined specifically.

4. Stochastic filtering. The system state one-step estimation is completed by the filtering algorithm and flowchart in section “Kalman filtering.”

5. Failure prognosis. If the system noise is assumed to be zero, and no new information is measured in equation (4), system state in future can be calculated by equation (5) in the probability sense.

Figure 4.

Failure prognosis process.

In this process, different uncertainties will be introduced and propagated step by step as shown in Figure 4. Next, the three types of uncertainties derived from different sources will be analyzed in this section.

Model parameter uncertainty

The real system cannot be modeled accurately, and the error of system mathematical model cannot be known in advance. For system state estimation and prediction based on stochastic filter, the system’s model is also completely unknown. Such error can be seen as system model parameter uncertainty. Taking the system description equation (1) as an example, system uncertainty cannot only be expressed by some parameters’ perturbation in system matrix ${\mathit{A}}_{\mathit{k}}$ but the whole matrix, which means that ${\mathit{A}}_{\mathit{k}}={\overline{\mathit{A}}}_{\mathit{k}}+\Delta {\mathit{A}}_{\mathit{k}}$ . The perturbation not only changes the distribution and number of zero-pole point of system but also changes the system structure. From viewpoint of theoretical analysis, such uncertainty should be described as unstructured set or structured set, which is suitable to be solved by some common analytical methods. ¹⁹

Measurement noise

Measurement noise or system output noise is the main uncertainty source during the failure prognosis process in a real environment. Environment noise, sensor inaccuracy, transmission noise, and data reduction effect can be seen as measurement noise. These noises can make the measurement result cannot reflect the true value of the measured item.

In terms of different mechanisms of the measurement noise, it can be classified as system error, random error, and gross error. System error is the measurement tool or system’s inherent error, and gross error is the error that is obviously inconsistent with the fact. Random error can be further divided as normal and abnormal distribution errors according to its PDF. In general case, system error and gross error can be removed by mathematical method, and the random error will be assumed to follow normal distribution. However, the statistical characteristic of dynamic noise in real environment is very hard to acquire exactly, and moreover, the initial “white” noise will also be colored as it is transferred through signal channel and preprocessed. This means, in most engineering applications, the measurement noise is only partly known or approximately or even totally unknown.^20,21

Failure threshold uncertainty

In the traditional overhaul or maintenance activities, the failure threshold, no matter such threshold value is given by human or intelligent monitor system, is given by experience and fixed in essence. However, it is just a traditional idea and may lead to the neglect of the uncertainty in engineering practice. Such assumption or neglect may not be unreasonable in most cases, as the “point” which the system performance degrades to failure, the failure threshold will be changed with the different operating environments, operators, and the system individual differences. ²² It is a challenging problem to consider the failure threshold’s uncertainty and model the real failure zone. In the sequel, the failure prognosis method with confidence interval will be presented in the next section, which considers the three kinds of uncertainties mentioned above.

Failure prognoses with multiple uncertainty sources

System parameter adaptive filtering

In general, the parameter uncertainty can be solved by the state augmentation; however, such method will increase the calculation, and sometimes, the estimation error will even become worse. For this reason, the idea of the separation of Friedland ²³ estimation was introduced in this section: the filtering process is separated robustly as the normal state estimation and system error estimation. For the state–space model in equation (1), it is assumed that there are some parameter perturbations in the system matrix which can be described as

$${\mathit{A}}_{\mathit{k}}={\overline{\mathit{A}}}_{\mathit{k}}+\Delta {\mathit{A}}_{\mathit{k}}$$

(7)

The uncertainty item $\Delta {\mathit{A}}_{\mathit{k}}$ satisfies

$$\Delta {\mathit{A}}_{\mathit{k}}={\mathit{E}}_{\mathit{k}}{\Delta }_{\mathit{k}}{\mathit{F}}_{\mathit{k}},\Vert {\Delta }_{\mathit{k}}\Vert \le 1$$

(8)

where ${\mathit{E}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{n}\times \mathit{r}}$ , ${\mathit{F}}_{\mathit{k}}\in {\mathbb{R}}^{\mathit{r}\times \mathit{n}}$ , and ${\Delta }_{\mathit{k}}\in {\mathbb{R}}^{\mathit{r}\times \mathit{r}}$ .

To simplify computation, the above three matrices are defined as

$${\Delta }_{\mathit{k}}=\mathit{diag}\left[{\delta }_1,{\delta }_2,\dots ,{\delta }_p\right],\mathrm{\hspace{0.17em}}{\mathit{E}}_{\mathit{k}}=\left[{\mathit{E}}_1,{\mathit{E}}_2,\dots ,{\mathit{E}}_p\right],\mathrm{\hspace{0.17em}}{\mathit{F}}_{\mathit{k}}={\left[{\mathit{F}}_1,{\mathit{F}}_2,\dots ,{\mathit{F}}_{\mathit{p}}\right]}^{\mathit{T}}$$

(9)

With equation (9), system model can be rewritten as

$${\mathit{x}}_{\mathit{k}}={\overline{\mathit{A}}}_{\mathit{k}-1}{\mathit{x}}_{\mathit{k}-1}+\left(\sum _{i=1}^p{\delta }_{\mathit{i}}{\mathit{E}}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}\right){\mathit{x}}_{\mathit{k}-1}+{\mathit{B}}_{\mathit{k}-1}{\mathit{v}}_{\mathit{k}-1}$$

(10)

Substituting equation (10) into measurement output expression in equation (1), the output error is

$${\hat{\mathit{w}}}_{\mathit{k}}={\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}={\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}[{\overline{\mathit{A}}}_{\mathit{k}-1}+\left(\sum _{i=1}^p{\delta }_{\mathit{i}}{\mathit{E}}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}\right)]{\hat{\mathit{x}}}_{\mathit{k}-1}$$

(11)

If the uncertainty item ${\delta }_{\mathit{i}}$ in equation (9) has its estimation value, this means ${\delta }_{\mathit{i}}={\hat{\delta }}_{\mathit{i}}+\Delta {\delta }_{\mathit{i}}$ , then equation (11) can be rewritten as

$${\hat{\mathit{w}}}_{\mathit{k}}={\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}[{\overline{\mathit{A}}}_{\mathit{k}-1}+\sum _{i=1}^p({\hat{\delta }}_{\mathit{i}}+\Delta {\delta }_{\mathit{i}}){\mathit{E}}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}]{\hat{\mathit{x}}}_{\mathit{k}-1}$$

(12)

Differentiate equation (12) with every uncertainty item ${\delta }_{\mathit{i}}$ to obtain

$$\frac{\partial {\hat{\mathit{w}}}_{\mathit{k}}}{\partial {\delta }_{\mathit{i}}}=-{\mathit{C}}_{\mathit{k}}{\mathit{E}}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}{\hat{\mathit{x}}}_{\mathit{k}-1},\mathit{i}=1,2,\dots ,\mathit{p}$$

(13)

$-{\mathit{C}}_{\mathit{k}}{\mathit{E}}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}$ can be rewritten as ${\mathit{H}}_{\mathit{k}}^{\delta }$ , that is, $\partial {\hat{\mathit{w}}}_{\mathit{k}}/\partial {\delta }_{\mathit{i}}={\mathit{H}}_{\mathit{k}}^{\delta }{\hat{\mathit{x}}}_{\mathit{k}-1}$ .

Set all the uncertainty items as

$${\delta }_{\mathit{k}}={\left[\begin{array}{cccc}\hfill {\delta }_1\hfill & \hfill {\delta }_2\hfill & \hfill \cdots \hfill & \hfill {\delta }_p\hfill \end{array}\right]}^{\mathit{T}}$$

Based on classical least square recursion, the system parameter update algorithm can be obtained as follows

$${\mathit{K}}_{\mathit{k}}^{\delta }={\mathit{P}}_{\mathit{k}}^{\delta }\mathit{H}{}_{\mathit{k}}^{\delta \mathit{T}}{({\mathit{R}}^{\delta }+{\mathit{H}}_{\mathit{k}}^{\delta }{\mathit{P}}_{\mathit{k}}^{\delta }\mathit{H}{}_{\mathit{k}}^{\delta \mathit{T}})}^{-1}$$

(14)

$${\delta }_{\mathit{k}+1}={\delta }_{\mathit{k}}+{\mathit{K}}_{\mathit{k}}^{\delta }{\hat{\mathit{w}}}_{\mathit{k}}$$

(15)

$${\mathit{P}}_{\mathit{k}+1}^{\delta }=(\mathit{I}-{\mathit{K}}_{\mathit{k}}^{\delta }{\mathit{H}}_{\mathit{k}}^{\delta }){\mathit{P}}_{\mathit{k}}^{\delta }$$

(16)

Based on measurement output error ${\hat{\mathit{w}}}_{\mathit{k}}$ in k th step, the deviation estimation process presented in equations (14)–(16) can be used to update system matrix’s uncertainty item $\Delta {\mathit{A}}_{\mathit{k}}={\mathit{E}}_{\mathit{k}}{\Delta }_{\mathit{k}}{\mathit{F}}_{\mathit{k}}$ with ${\delta }_{\mathit{k}}$ , which realizes system parameter’s adaptive updating.

Robust filters design for unknown measurement noise

System or measurement noise in Kalman filtering process is defined as variance matrices ${\mathit{Q}}_{\mathit{k}}$ and ${\mathit{R}}_{\mathit{k}}$ in equation (2). In Figure 1, such prior information is assumed to be known and used to compute filtering gain ${\mathit{K}}_{\mathit{k}}$ and estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ . The filtering gain ${\mathit{K}}_{\mathit{k}}$ decides the weight between new information ${\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}$ and old information ${\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}$ in state estimation equation as follows: ${\hat{\mathit{x}}}_{\mathit{k}}={\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}({\mathit{y}}_{\mathit{k}}-{\mathit{C}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1})$ , and the gain ${\mathit{K}}_{\mathit{k}}$ is used to compute the estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ in equation (4). Furthermore, the estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ decides the estimation error upper bound in every calculation step; meanwhile, such item also has relationship with one-step forecasting mean-squared error ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ .

In all, the gain calculation loop, as shown in Figure 1, which contains the three items (i.e. ${\mathit{K}}_{\mathit{k}}$ , ${\mathit{P}}_{\mathit{k}}$ , and ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ ), is the core of Kalman filtering process. So, if we can find some gain ${\mathit{K}}_{\mathit{k}}$ which guarantees that the filter’s one-step forecasting mean-squared error ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ and estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ are bounded all the time, or in other words, state estimation error in all calculation steps is acceptable, and the filter is robust. For the system in equation (1), its one-step forecasting mean-squared error ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ is

$${\mathit{P}}_{\mathit{k}/\mathit{k}-1}={\mathit{A}}_{\mathit{k}-1}{\mathit{P}}_{\mathit{k}-1}{\mathit{A}}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{B}}_{\mathit{k}-1}{\mathit{QB}}_{\mathit{k}-1}^{\mathit{T}}$$

Substituting the estimation mean-squared error ${\mathit{P}}_{\mathit{k}-1}$ in former calculation step and one-step forecasting mean-squared error ${\mathit{P}}_{\mathit{k}-1/\mathit{k}-2}$ in equation (4) in the above equation and repeating this process, we can obtain ²⁴

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{P}}_{\mathit{k}/\mathit{k}-1}}& =\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{1}{\Pi }}{\Lambda }_{\mathit{i}}\right){\mathit{A}}_0{\mathit{P}}_0{\mathit{A}}_0^{\mathit{T}}{\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{1}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}\hfill \\ \multicolumn{1}{c}{}& +\sum _{\mathit{j}=1}^{\mathit{k}-1}\left[\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right){\mathit{B}}_{\mathit{j}-1}\mathit{Q}{\mathit{B}}_{\mathit{j}-1}^{\mathit{T}}{\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}\right]\hfill \\ \multicolumn{1}{c}{}& +\sum _{\mathit{j}=\mathit{k}-1}^2\left[\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right){\mathit{A}}_{\mathit{j}-1}{\mathit{K}}_{\mathit{j}-1}\mathit{R}{\mathit{K}}_{\mathit{j}-1}^{\mathit{T}}{\mathit{A}}_{\mathit{j}-1}^{\mathit{T}}{\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}\right]\hfill \\ \multicolumn{1}{c}{}& +{\mathit{A}}_{\mathit{k}-1}{\mathit{K}}_{\mathit{k}-1}\mathit{R}{\mathit{K}}_{\mathit{k}-1}^{\mathit{T}}{\mathit{A}}_{\mathit{k}-1}^{\mathit{T}}+{\mathit{B}}_{\mathit{k}-1}\mathit{Q}{\mathit{B}}_{\mathit{k}-1}^{\mathit{T}}\hfill \end{array}$$

(17)

where ${\Lambda }_{\mathit{i}}={\mathit{A}}_{\mathit{i}}(\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})$ . Then, the one-step forecasting mean-squared error’s upper bound can be obtained as

$$\begin{array}{cc}\multicolumn{1}{c}{\Vert {\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\Vert }_2=}& \Vert \left(\underset{\mathit{i}=\mathit{k}-1}{\overset{1}{\Pi }}{\Lambda }_{\mathit{i}}\right){\Vert }_2\Vert {\mathit{A}}_0{\mathit{P}}_0{\mathit{A}}_0^{\mathit{T}}{\Vert }_2\Vert {\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{1}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}{\Vert }_2+\sum _{\mathit{j}=1}^{\mathit{k}-1}\left[{\Vert \left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)\Vert }_2{\Vert {\mathit{B}}_{\mathit{j}-1}\mathit{Q}{\mathit{B}}_{\mathit{j}-1}^{\mathit{T}}\Vert }_2{\Vert {\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}\Vert }_2\right]\hfill \\ \multicolumn{1}{c}{}& +\sum _{\mathit{j}=\mathit{k}-1}^2\left[{\Vert \left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)\Vert }_2{\Vert {\mathit{A}}_{\mathit{j}-1}{\mathit{K}}_{\mathit{j}-1}\Vert }_2{\Vert \mathit{R}\Vert }_2{\Vert {\mathit{K}}_{\mathit{j}-1}^{\mathit{T}}{\mathit{A}}_{\mathit{j}-1}^{\mathit{T}}\Vert }_2{\Vert {\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}{\Lambda }_{\mathit{i}}\right)}^{\mathit{T}}\Vert }_2\right]\hfill \\ \multicolumn{1}{c}{}& +\Vert {\mathit{A}}_{\mathit{k}-1}{\mathit{K}}_{\mathit{k}-1}{\Vert }_2\Vert \mathit{R}{\Vert }_2\Vert {\mathit{K}}_{\mathit{k}-1}^{\mathit{T}}{\mathit{A}}_{\mathit{k}-1}^{\mathit{T}}{\Vert }_2+\Vert {\mathit{B}}_{\mathit{k}-1}{\mathit{QB}}_{\mathit{k}-1}^{\mathit{T}}{\Vert }_2\hfill \end{array}$$

(18)

Let ${\lambda }_1^{{\Lambda }_{\mathit{i}}}$ denote the maximum eigenvalue of ${\Lambda }_{\mathit{i}}^{\mathit{T}}{\Lambda }_{\mathit{i}}$ and ${\lambda }_1^{{\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}}$ denote the maximum eigenvalue of $({\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}{)}^{\mathit{T}}{\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}$ , then based on the two-norm definition $\Vert {\Lambda }_{\mathit{i}}{\Vert }_2=\sqrt{{\lambda }_1^{{\Lambda }_{\mathit{i}}}}$

$$\begin{array}{c}\Vert {\mathit{P}}_{\mathit{k}/\mathit{k}-1}{\Vert }_2\le \underset{\mathit{i}=\mathit{k}-1}{\overset{1}{\Pi }}\left|{\lambda }_1^{{\Lambda }_{\mathit{i}}}\right|\Vert {\mathit{A}}_0{\mathit{P}}_0{\mathit{A}}_0^{\mathit{T}}{\Vert }_2\hfill \\ +\sum _{\mathit{j}=1}^{\mathit{k}-1}\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}\left|{\lambda }_1^{{\Lambda }_{\mathit{i}}}\right|{\Vert {\mathit{B}}_{\mathit{j}-1}\mathit{Q}{\mathit{B}}_{\mathit{j}-1}^{\mathit{T}}\Vert }_2\hfill \\ +\sum _{\mathit{j}=\mathit{k}-1}^2\left\{\left(\underset{\mathit{i}=\mathit{k}-1}{\overset{\mathit{j}}{\Pi }}\left|{\lambda }_1^{{\Lambda }_{\mathit{i}}}\right|\right){\lambda }_1^{{\mathit{A}}_{\mathit{j}-1}{\mathit{K}}_{\mathit{j}-1}}\right\}\Vert \mathit{R}{\Vert }_2\hfill \\ +\left|{\lambda }_1^{{\mathit{A}}_{\mathit{k}-1}{\mathit{K}}_{\mathit{k}-1}}\right|\Vert \mathit{R}{\Vert }_2+\Vert {\mathit{B}}_{\mathit{k}-1}Q{\mathit{B}}_{\mathit{k}-1}^{\mathit{T}}{\Vert }_2\hfill \end{array}$$

(19)

If the spectral radius $\left|{\lambda }_1^{{\Lambda }_{\mathit{i}}}\right|=\rho ({\Lambda }_{\mathit{i}}^{\mathit{T}}{\Lambda }_{\mathit{i}})<1$ and, meanwhile, $\rho [({\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}{)}^{\mathit{T}}{\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}]\le \alpha$ , then equation (19) is bounded. With the estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ ’s expression in equation (4), if the result of one-step forecasting mean-squared error ${\mathit{P}}_{\mathit{k}/\mathit{k}-1}$ is substituted into its expression, then the similar conclusion for ${\mathit{P}}_{\mathit{k}}$ can be obtained as follows: If the spectral radius $|{\lambda }_1^{(1-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})}|=\rho [(1-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}{)}^{\mathit{T}}(1-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})]<\beta$ and, meanwhile, $\rho [({\mathit{K}}_{\mathit{k}}{)}^{\mathit{T}}{\mathit{K}}_{\mathit{k}}]\le \gamma$ , then estimation mean-squared error ${\mathit{P}}_{\mathit{k}}$ is bounded. To obtain the numerical resolution of the four conditions mentioned above, the robust filtering theorem is presented in the following.

Theorem 1

The filter with filtering gain $\mathit{K}$ is robust if there exist real constants $\alpha ,\beta ,\mathrm{and}\gamma$ and the gain which satisfy the following linear matrix inequality (LMIs)

$$\left[\begin{array}{cc}\hfill {{\mathit{A}}_{\mathit{i}}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}-{{\mathit{A}}_{\mathit{i}}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{{\mathit{K}}_{\mathit{i}}}^{\mathit{T}}{{\mathit{A}}_{\mathit{i}}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}-\mathit{I}\hfill & \hfill {\mathit{C}}_{\mathit{i}}^{\mathit{T}}{{\mathit{K}}_{\mathit{i}}}^{\mathit{T}}{{\mathit{A}}_{\mathit{i}}}^{\mathit{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}}& \hfill -\mathit{I}\hfill \end{array}\right]<0$$

(20)

$$\left[\begin{array}{cc}\hfill -{\alpha }^2\mathit{I}\hfill & \hfill {{\mathit{K}}_{\mathit{i}}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}^{\mathit{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}}& \hfill -\mathit{I}\hfill \end{array}\right]\le 0$$

(21)

$$\left[\begin{array}{cc}\hfill \mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}-{\beta }^2\mathit{I}\hfill & \hfill {\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}}& \hfill -\mathit{I}\hfill \end{array}\right]\le 0$$

(22)

$$\left[\begin{array}{cc}\hfill -{\gamma }^2\mathit{I}\hfill & \hfill {{\mathit{K}}_{\mathit{i}}}^{\mathit{T}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathit{i}}}& \hfill -\mathit{I}\hfill \end{array}\right]\le 0$$

(23)

Proof

With the above derivation, we can conclude that the state estimation error in all steps is acceptable if the following four conditions are valid

$$\rho \left({\Lambda }_{\mathit{i}}^{\mathit{T}}{\Lambda }_{\mathit{i}}\right)<1$$

$$\rho \left[{\left({\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}\right)}^{\mathit{T}}{\mathit{A}}_{\mathit{k}}{\mathit{K}}_{\mathit{k}}\right]\le \alpha$$

$$\rho \left[{(1-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})}^{\mathit{T}}(1-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})\right]<\beta$$

$$\rho \left[{\left({\mathit{K}}_{\mathit{k}}\right)}^{\mathit{T}}{\mathit{K}}_{\mathit{k}}\right]\le \gamma$$

These conditions can be transferred to

$${\left[{\mathit{A}}_{\mathit{i}}(\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})\right]}^{\mathit{T}}\left[{\mathit{A}}_{\mathit{i}}(\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})\right]<\mathit{I}$$

$${\left({\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}\right)}^{\mathit{T}}\left({\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}\right)\le \alpha \mathit{I}$$

$${(\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})}^{\mathit{T}}(\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}})\le \beta \mathit{I}$$

$${\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}\le \mathit{I}$$

and they are equal to

$${\mathit{A}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}-{\mathit{A}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}+{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-\mathit{I}<0$$

$${\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}^{\mathit{T}}{\mathit{A}}_{\mathit{i}}{\mathit{K}}_{\mathit{i}}-\alpha \mathit{I}\le 0$$

$$\mathit{I}-{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}+{\mathit{C}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}-\beta \mathit{I}\le 0$$

$${\mathit{K}}_{\mathit{i}}^{\mathit{T}}{\mathit{K}}_{\mathit{i}}-\gamma \mathit{I}\le 0$$

With the Shur formula, ²⁰ LMIs (20)–(23) can be obtained.

Remarks 1

According to our analysis, as the real constants $\alpha ,\beta ,\mathrm{and}\gamma$ become greater, the convergence speed of filter will become slower or even divergent. However, if such parameters become too small, LMIs will become unsolvable.

Failure prognosis with failure threshold distribution

Based on the advanced robust filtering on section “Robust filters design for unknown measurement noise,” the forward failure prognosis value calculation can be started from “current filtering state” ${\hat{\mathit{x}}}_{{\mathit{k}}_0}$ . Similar to equation (5), step-by-step extrapolation process can be rewritten as

$${\hat{\mathit{x}}}_{{\mathit{k}}_{\mathit{i}+1}}={\overline{\mathit{A}}}_{{\mathit{k}}_0}{\hat{\mathit{x}}}_{{\mathit{k}}_{\mathit{i}}}$$

(24)

$${\mathbf{P}}_{{\mathit{k}}_{\mathit{i}+1}}={\overline{\mathit{A}}}_{{\mathit{k}}_{\mathit{i}}}{\mathbf{P}}_{{\mathit{k}}_{\mathit{i}}}{\overline{\mathit{A}}}_{{\mathit{k}}_{\mathit{i}}}^{\mathit{T}}\mathit{i}=0,\dots ,\mathit{n}$$

(25)

Since there is no variance item in robust filter in section “Robust filters design for unknown measurement noise,” the variance matrix ${\mathit{P}}_{{\mathit{k}}_0}$ in “current time state” should be reconstructed as

$${\hat{\mathit{x}}}_{{\mathit{k}}_0}={\hat{\mathit{x}}}_{{\mathit{k}}_0/{\mathit{k}}_0-1}+{\mathit{K}}_{{\mathit{k}}_0}({\mathit{y}}_{{\mathit{k}}_0}-{\hat{y}}_{{\mathit{k}}_0})$$

$${\mathit{P}}_{{\mathit{k}}_0}=\mathit{E}\left[{\tilde{\mathit{x}}}_{{\mathit{k}}_0}{\tilde{\mathit{x}}}_{{\mathit{k}}_0}^{\mathit{T}}\right]=\left[{\mathit{K}}_{{\mathit{k}}_0}({\mathit{y}}_{{\mathit{k}}_0}-{\hat{\mathit{y}}}_{{\mathit{k}}_0})\right]{\left[{\mathit{K}}_{{\mathit{k}}_0}({\mathit{y}}_{{\mathit{k}}_0}-{\hat{\mathit{y}}}_{{\mathit{k}}_0})\right]}^{\mathit{T}}$$

As shown in Figure 3 in section “Failure prognosis approach,” the cumulative probability of failure before time $({\mathit{k}}_{\mathit{i}}-{\mathit{k}}_0)\mathit{T}$ can be integrated by equation (6) as

$$\mathrm{Pr}({\mathit{T}}_{\mathit{f}}\le ({\mathit{k}}_{\mathit{i}}-{\mathit{k}}_0)\mathit{T})=1-\Phi \left(\frac{{\mathit{x}}_{\mathit{th}}-{\hat{\mathit{x}}}_{{\mathit{k}}_{\mathit{i}}}}{\sqrt{{\mathbf{P}}_{{\mathit{k}}_{\mathit{i}}}(1,1)}}\right)$$

(26)

where $\mathit{T}$ is the sampled time interval, ${\mathit{x}}_{\mathit{th}}$ is the predefined failure threshold, and the failure prognosis value’s PDF can be obtained by the above equation’s derivation.

If the confidence probability is set to $\zeta \%$ , then

$$1-\Phi \left(\frac{{\mathit{x}}_{\mathit{th}}-{\widehat{\mathit{x}}}_{{\mathit{k}}_{\mathit{i}}}}{\sqrt{P_{{\mathit{k}}_{\mathit{i}}}\left(1,1\right)}}\right)=\zeta \%\Rightarrow \Phi \left(\frac{{\mathit{x}}_{\mathit{th}}-{\widehat{\mathit{x}}}_{{\mathit{k}}_{\mathit{i}}}}{\sqrt{P_{{\mathit{k}}_{\mathit{i}}}\left(1,1\right)}}\right)=1-\zeta \%$$

(27)

With equations (24), (25), and (27), the failure prognosis value with confidence probability $\zeta \%$ can be obtained as ${\mathit{T}}_{\mathit{f}}=({\mathit{k}}_{\mathit{i}}-{\mathit{k}}_0)\mathit{T}$ .

In general, the actual failure threshold ${\mathit{x}}_{\mathit{th}}$ will vary with the different environments. Based on this, the failure threshold ${\mathit{x}}_{\mathit{th}}$ is defined to follow the probability distribution function $\mathit{f}({\mathit{x}}_{\mathit{th}})$ in this article and then equation (26) can be rewritten by total probability formula

$$\mathrm{Pr}\left({\mathit{T}}_f\le \left({\mathit{k}}_i-{\mathit{k}}_0\right)\mathit{T}\right)={\displaystyle \underset{{\mathit{x}}_{\mathit{th}}^{low}}{\overset{x_{th}^{upp}}{\int }}\left\{f\left(x_{th}\right)\left[1-\Phi \left(\frac{{\mathit{x}}_{\mathit{th}}-{\widehat{\mathit{x}}}_{{\mathit{k}}_i}}{\sqrt{P_{k_i}\left(1,1\right)}}\right)\right]\right\}}\mathit{d}{\mathit{x}}_{\mathit{th}}$$

(28)

${\mathit{x}}_{\mathit{th}}^{\mathit{upp}}$ and ${\mathit{x}}_{\mathit{th}}^{\mathit{low}}$ are the upper and lower bounds of function $\mathit{f}({\mathit{x}}_{\mathit{th}})$ .

Remark 2

If the failure threshold follows normal distribution ${\mathit{x}}_{\mathit{th}}\in \mathit{N}({\mu }_{\mathit{th}},{\sigma }_{\mathit{th}}^2)$ , then we can choose the upper and lower bounds of threshold as ${\mu }_{\mathit{th}}\pm 3{\sigma }_{\mathit{th}}$ , which guarantee the interval confidence up to 99.73%.

Failure prognosis value estimation algorithms

In section “Failure prognoses with multiple uncertainty sources,” the probability of system state falling into the failure threshold $({\mathit{x}}_{\mathit{th}}^{\mathit{upp}},{\mathit{x}}_{\mathit{th}}^{\mathit{low}})$ is given by equation (28). This means, if we can measure the key parameters and estimate the future state’s value, then the failure prognosis can be carried out to help engineers or managers about spare parts management or maintenance plans development. Based on the failure prognosis equations in section “Failure prognoses with multiple uncertainty sources,” the failure prognostic algorithm can be summarized and listed as follows:

• Step 1. Initialize the system state ${\mathit{x}}_0$ , system matrix variance matrix ${\mathit{P}}_0^{\delta }$ , and matrix noise matrix ${\mathit{R}}^{\delta }$ .

• Step 2. Calculate one-step forecasting ${\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}={\overline{\mathit{A}}}_{\mathit{k}-1}{\hat{\mathit{x}}}_{\mathit{k}-1}$ and output estimation ${\hat{\mathit{y}}}_{\mathit{k}}={\overline{\mathit{C}}}_{\mathit{k}}{\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}$ based on it.

• Step 3. Calculate the filter gain ${\mathit{K}}_{\mathit{k}}$ in current state based on the LMIs in theorem 1.

• Step 4. Get the state estimation as ${\hat{\mathit{x}}}_{\mathit{k}}={\hat{\mathit{x}}}_{\mathit{k}/\mathit{k}-1}+{\mathit{K}}_{\mathit{k}}({\mathit{y}}_{\mathit{k}}-{\hat{\mathit{y}}}_{\mathit{k}})$ .

• Step 5. Assume the initial system matrix parameter error is 0 and then update parameter by the parameter identification algorithm in equations (13)–(16).

• Step 6. Reset ${\overline{\mathit{A}}}_{\mathit{k}}={\overline{\mathit{A}}}_{\mathit{k}-1}+\delta {\mathit{A}}_{\mathit{k}-1}$ .

• Step 7. Predict the system’s future state value step by step by equations (24) and (25) and then complete the failure prognosis with predefined failure threshold distribution and failure probability.

Such process can be shown in Figure 5.

Figure 5.

Failure prognosis value estimation procedure.

Simulation examples

In this section, the advanced failure prognosis method proposed in this article will be applied to an aircraft fuel feed system simulation case. First, the fuel feed system’s failure modes are analyzed and then the proposed method and other classical methods are used to complete failure prognosis and compare with each other.

Failure effect analysis of fuel feed system

The function of fuel system is to continuously and stably provide engine and auxiliary power unit (APU) with adequate fuel under permitted working conditions and flight states. Modern aircraft fuel feed system consists of electronic control unit, fuel pipeline, and different valves and is a hybrid system with mechanical, electronic, and liquid component, which is shown in Figure 6.

Figure 6.

Aircraft fuel feed system diagram.

Each fuel tank has two booster pumps; fuel booster pump extracts fuel and supplies it to the engine feed manifold by pipeline. According to the fuel system operating principle, the component’s performance or reliability on supply line is important to guarantee fuel supply. Booster pump is the key component on supply line, which provides fuel press in most flight stages. In this article, we take the booster pump subsystem as an example.

Booster pump boosts the fuel by its vane’s high-speed rotation, which will introduce cavitation erosion. Under the high-speed flowing and pressure variation condition, the air bubble will be introduced on the vane’s surface as the pressure on the vane is lower than fluid’s steam pressure. And meanwhile, the gas dissolved in fluid may be also separated out to air bubble. After this, the air bubble will be broken as the fluid’s pressure is greater than bubble’s pressure; such sudden action will make huge temperature impact on vane’s surface, which leads to the vane’s fatigue and breaking off or even big pit. Such phenomenon is called cavitation erosion, as shown in Figure 7. Cavitation erosion will decrease the vane area and then makes the pump’s driving force weaker and weaker. Meanwhile, the output pressure will also become smaller and smaller, until it cannot satisfy fuel feed system’s performance index. This means, fuel supply pressure is sensitive to the inner failure mode and can be used to monitor system’s health state.

Figure 7.

Vane’s cavitation erosion.

Based on the characteristic parameter simulation output, we can complete the fuel feed system’s failure prognosis with the proposed method. Figure 8 shows the output pressure value of hydraulic booster pump in aircraft fuel feed system with different wear parameters during the whole life period; the x-axis represents system running time, and y-axis shows the pressure values. As shown in Figure 8, during the whole operation cycle, pump’s output pressure value variation presents different degradation trends with different wear conditions. In this figure, the x-axis representing time is operational time (×10⁶ s).

Figure 8.

Pressure curves during the whole life period under different wear coefficient conditions.

Parameter adaptive filtering results

As shown in Figure 8, fuel system’s output pressure varies dramatically with the different wear levels. This means, system mathematical model which describes the degradation process also has to be changed; in other words, the actual system’s model has uncertainty factors.

From viewpoint of theoretical model, such variation can be described as perturbation $\Delta {\mathit{A}}_{\mathit{k}}$ in equation (1). Then, the system parameter adaptive filtering method in section “Failure prognoses with multiple uncertainty sources” can be employed to estimate system state in the following. In order to describe the system parameter uncertainty, the system’s motion is assumed to satisfy classical constant acceleration model ²⁵

$$\stackrel{·}{\mathit{a}}=\eta$$

(29)

where $\mathit{a}$ is the system motion acceleration. For uncertainty description, system input $\eta$ can be written as

$$\eta =\mathit{b}+{\eta }_{\mathit{w}}$$

where $\mathit{b}$ is the uncertainty item in system matrix and ${\eta }_{\mathit{w}}$ is the system’s random noise. All this error or uncertainty will be changed as environment or workload varies and cannot be known beforehand, which forms the system model’s uncertainty $\Delta {\mathit{A}}_{\mathit{k}}$ .

In contrast, the classical Kalman filer or nonparameter adaptive filter and parameter adaptive filter are adopted to estimate the actual output pressure data with different wear levels. Figures 9 and 10 and Table 1 illustrate the estimation results’ comparison with different methods.

Figure 9.

Comparison of estimation error with different wear levels.

Figure 10.

Comparison of estimation results with given wear level.

Table 1.

Mean-squared error of two estimation methods with different wear levels.

	Adaptive filter	Nonadaptive filter
Wear level 1	0.000168453	0.05604427
Wear level 2	0.000224	0.065535
Wear level 3	0.016783	0.070623
Wear level 4	0.011715	0.066545
Wear level 5	0.01424	0.214038
Wear level 6	0.015139	0.114475
Wear level 7	0.018104	0.1254

Figures 9 and 10 show that the estimation error of nonadaptive filter obviously becomes greater as system workload increases or wear level becomes worse, but the adaptive filter which updates the system parameter in each estimation step is fit for different wear levels and can reduce the estimation error. Obviously, the subsequent failure prognosis case simulation in section “Failure prognosis estimation results” will be affected by these different estimation effects dramatically.

Robust filtering results

The model in equation (29) with unknown noise characteristic or colored noise will be used to verify robust filtering method’s performance. As discussed at the end of section “Failure prognoses with multiple uncertainty sources,” the three parameters $\alpha ,\beta ,\mathrm{and}\gamma$ in equations (21)–(23) decide the filter’s convergence speed and the estimation mean-squared error matrix’s upper bound. In this section, we choose $\alpha ,\beta ,\mathrm{and}\gamma$ as 1.2; then, the estimation curves can be obtained as shown in Figures 11 and 12 for colored noise and unknown noise characteristic, respectively.

Figure 11.

Comparison of absolute estimation error with colored noise condition.

Figure 12.

Comparison of absolute estimation error with unknown noise condition.

Figure 11 shows the absolute estimation error of the proposed robust filter and classical Kalman filter when system output pressure data have colored noise. With unknown characteristic noise in measured data, Figure 12 illustrates the absolute estimation error for four different filters, which are classical Kalman filter with known noise level, robust filter, classical Kalman filter with unknown noise level, and the H_∞ filter. From Figures 11 and 12, the following is clear that

1. The robust filter is insensitive to colored noise and can guarantee system estimation error better than classical Kalman filter.

2. The robust filter’s performance is same as classical Kalman filter with known noise level and may be a little better than H_∞ filter.

3. Obviously, robust filter is better than the Kalman filter with unknown noise level.

In summary, for system with unknown noise characteristic, such robust approach can improve estimation efficiency. And the theoretical analysis in section “Failure prognoses with multiple uncertainty sources” is verified by the simulation results.

Failure prognosis estimation results

The approach derived in section “Failure prognoses with multiple uncertainty sources” will be used for failure prognosis of fuel booster pump system. Through the simulation in section “Parameter adaptive filtering results,” we can conclude that the model exhibits good tracking ability for all the testing data, and one of the results is shown in Figure 10.

The adaptive and robust characteristics of the above approach can be used to update the system parameters in the presence of new information and are insensitive to unknown noise. The failure prognosis results with equations (24)–(27) and the deterministic threshold at several time points are plotted in Figures 13–15. In our simulation, the deterministic threshold for booster pump is set at 2 bar, and the confidence probability is 50%.

Figure 13.

Comparison of advanced and classical failure prognosis results with single threshold.

Figure 14.

Comparison of advanced and classical failure prognosis PDF results at different time point 1.

Figure 15.

Comparison of advanced and classical failure prognosis PDF results at different time point 2.

In Figure 13, the solid line is the measured output signal of booster pump pressure, the dash line is real decay parameter’s value, the dot line is parameter prediction zone by advanced failure prognosis estimation method in this article, and the dash–dot line is parameter prediction zone by classical failure prognosis estimation method which is given in section “Fault prognosis theory based on stochastic filtering.” In different time points, the prediction result of our method is fit better than the classical method, and moreover, as shown in Figure 13, the prediction zone of classical method has deviated from the actual failure time point.

Figure 14 illustrates the conditional PDFs for failure prognosis obtained at each time point with different methods, respectively; the bold solid line is the prediction results for advanced method in this article and dash line is for classical method, and meanwhile, the dots represent the actual failure prognosis at the corresponding time point. The two PDFs are obviously different, and the solid ones are narrower than the dash ones which are obtained from the classical method.

Figure 15 shows the plane information about the failure prognosis in different time points. The x-axis indicates the booster pump’s working time, and the y-axis indicates the system failure time. The dot line throughout the upper left corner to the lower right corner in Figure 15 represents the actual system failure time in current time, and the thick solid line is obtained by the failure prognosis points set segment connection, whose points were obtained by the proposed method in this article; similarly, the thick dashed line’s points were obtained by the classical method. In order to show clearly, the PDFs of each failure prognosis point with different methods are given by the dashed line. Obviously, the curves show that the proposed method in this article is more reliable, and the solid line is more fit to the actual failure prognosis line.

In summary, with parameter updating and robust treatment of system estimation in each step, the proposed failure prognosis method which considers multiple uncertainty sources shows better failure prognosis efficiency and accuracy.

From the analysis in section “Uncertainty source analysis,” the actual component’s failure time point always belongs to a failure time interval in the sense of probability due to the complex field environment or device individual difference, rather than a single threshold just like the above simulation, which means that the actual failure threshold may follow some function distribution $\mathit{f}({\mathit{x}}_{\mathit{th}})$ . In this section, the failure threshold is assumed to follow a normal distribution around the single threshold 2 bar, and the variance is 0.1, for example, $\mathit{N}(2,0.1)$ . And then the failure prognosis PDF calculation by both the single threshold and threshold interval with distribution at different time points can be plotted in Figure 16.

Figure 16.

Comparison of failure prognosis PDF results with single threshold and threshold zone.

Similar to Figure 14, the dots represent the actual failure prognosis at the corresponding time point, and of course, the failure prognoses by the two methods which consider or do not consider threshold distribution contain such points. As current time point moving, the curves become narrow and the failure prognosis results become more credible. For the first curve set, the two results near the time 40 (bold solid line for single threshold and dash line for another) are almost the same because of failure time is far from current time, and there is only a little information available for prediction. However, as the time point moving, the two results become more different gradually. Figure 16 shows that compared with the failure prognosis with single threshold, the introduction of threshold distribution will reduce component’s unnecessary early replacement.

Conclusion

In this article, a failure prognosis approach with multiple uncertainty sources has been presented. Such approach utilizes adaptive updating scheme, robust filter, and total probability formula to form the advanced failure prognosis tool which can predict the residual life based on the measured data. The fault prognostic process was listed in section “Failure prognosis value estimation algorithms” to simplify the actual use of such approach.

Simulated results for the fuel feed system’s booster pump showed the validity of the proposed method. Further research should be focused on future usage uncertainty which has relationship with maintenance decision, and another meaningful area is the real engineering environment, which will perform its significant value in maintenance or logistics support.Furthermore, we would like to see the proposed failure prognosis approach can be used in flight control problems.^26–28