Application of the Kalman Filter in Monitoring, Diagnosis, and Fault Parrying Problems for Observable Dynamical Systems

Abstract

This section is devoted to approaches for the formation, estimation and application of diagnostic parameters, as well as stochastic models in monitoring problems. The procedures for increasing the reliability and depth of dynamic systems monitoring are substantiated based on the decomposition of diagnostic models, the use of combined statistical criteria and the processing of observations in “forward” and “backward” time using the extended Kalman filter. The possibility of detecting and parrying anomalous observations using robust filtering procedures is shown. Such procedures are based on the use of an influence function that establishes the level of confidence for formed observations. U-D modification of the Kalman filter with an influence function in the observation selection loop is presented. Here U is an upper triangular matrix with unit diagonal elements and D is a diagonal matrix. Results of the mathematical simulation are given.

1. Introduction

At present, the problem of increasing the reliability of dynamical systems (DSs), which are part of objects of various physical natures, still remains topical. The solution to the above problem can be based on the application of analytical approaches to fault detection and isolation (FDI). In observable dynamic systems, the FDI problem is solved using external information, diagnostic models, and decision rules. When constructing diagnostic models, a mathematical description of the reference (unperturbed) DS and the actual (perturbed) DS is used. The reference vector Y(t) and the actual vector $Y_a\left(t\right)$ of state parameters correspond to such DSs. The dynamics of changes in these vectors are described by the following differential equations:

$$\mathit{dY}\left(t\right)/\mathit{dt}=\dot{Y}\left(t\right)=F\left[Y\left(t\right)\right],$$E1

$${\dot{Y}}_a\left(t\right)=F\left[Y_a\left(t\right)\right]+G\left(t\right)\mathrm{\xi }\left(t\right),$$E2

where F(…) is a matrix of variable coefficients characterizing the dynamics of changes in DS parameters; $\mathrm{\xi }\left(t\right)$ is a vector of random deviations of DS parameters from the required values; $E\left[\mathrm{\xi }\left(t\right){\mathrm{\xi }}^T\left(t-\mathrm{\tau }\right)\right]=Q\left(t\right)\mathrm{\delta }\left(t-\mathrm{\tau }\right)$ is the covariance matrix; Q(t) is a matrix of weighting coefficients; $\mathrm{\delta }\left(t-\mathrm{\tau }\right)$ is the delta-function; E […] is the operator of mathematical expectation; $G\left(t\right)$ is a matrix of variable coefficients characterizing the dynamics of changes in random deviations.

Taking into account relations (1) and (2), the equation for DS errors will have the form

$$\mathit{dx}\left(t\right)/\mathit{dt}=\dot{x}\left(t\right)=A\left(t\right)x\left(t\right)+G\left(t\right)\mathrm{\xi }\left(t\right),$$E3

where $x\left(t\right)=\mathrm{\Delta }Y\left(t\right)=Y_a\left(t\right)-Y\left(t\right)$ is the vector of DS errors; $A\left(t\right)={\left.\partial F\left[Y\left(t\right)\right]/\partial Y\right|}_{Y\left(t\right)=Y_a\left(t\right)}$ is a matrix of variable coefficients characterizing the dynamics of changes in DS errors.

The estimates $\widehat{x}\left(t\right)$ of DS errors can be obtained by processing observations z(t) using the extended Kalman filter (EKF) [1]:

$$z\left(t\right)=h\left[Y_a\left(t\right)\right]-h{\left[Y\left(t\right)\right]}_{\mathrm{SEI}},$$E4

where SEI is the sensor of external information; $h{\left[Y\left(t\right)\right]}_{\mathrm{SEI}}$ is the observable parameters formed by SEI and having a model $h{\left[Y\left(t\right)\right]}_{\mathrm{SEI}}=h\left[Y\left(t\right)\right]+\mathrm{\vartheta }\left(t\right)$; and $\mathrm{\vartheta }\left(t\right)$ is a vector of random deviations of SEI parameters from the reference with covariance matrix $E\left[\mathrm{\vartheta }\left(t\right){\mathrm{\vartheta }}^{\mathrm{T}}\left(t-\mathrm{\tau }\right)\right]=R\left(t\right)\mathrm{\delta }\left(t-\mathrm{\tau }\right)$.

The mathematical model of observation (4) has the form:

$$z\left(t\right)=H\left(t\right)x\left(t\right)+\mathrm{\vartheta }\left(t\right),$$E5

where $H\left(t\right)={\left.\partial h\left[Y\left(t\right)\right]/\partial Y\right|}_{Y\left(t\right)=Y_a\left(t\right)}$ is the matrix for the relation of observed parameters and the vector of DS errors.

The extension of the traditional linear Kalman filter is associated with its application in optimization problems of nonlinear dynamic systems of the form (2)–(5). The observable DS with EKF in the errors estimation loop is shown in Figure 1. $\widehat{Z}$ is the predicted estimate of observations; ^ is the symbol for estimate; and $\mathrm{\nu }$ is the vector of residuals. Such a vector can be used to form diagnostic parameters for monitoring observed dynamical systems with EKF in the fault detection loop.

Figure 1.

An observable dynamical system with the EKF in an error estimation loop.

Discrete observations $Z_i=Z\left(t_i\right)$ are related with the following DS error equation:

$$x_i={\mathrm{\Phi }}_i{x}_{i-1}+{\mathrm{\Gamma }}_i{\mathrm{\xi }}_{i-1},$$E6

where $x_i=x\left(t_i\right)$; ${\mathrm{\Phi }}_i$ matrix is the solution of the homogeneous part of Eq. (3).

$$d\mathrm{\Phi }\left(t,t_{i-1}\right)/\mathit{dt}={\dot{\mathrm{\Phi }}}_i=A\left(t\right)\mathrm{Ф}\left(t,t_{i-1}\right);\mathrm{Ф}\left(t_{i-1},t_{i-1}\right)=I;$$

I is an identity matrix; ${\mathrm{\Gamma }}_i$ matrix is the solution of the inhomogeneous part of Eq. (3):

$$d\mathrm{\Gamma }\left(t,t_i\right)/\mathit{dt}={\dot{\mathrm{Г}}}_i=A\left(t\right)\mathrm{Г}\left(t,t_{i-1}\right)+G\left(t\right);\mathrm{\Gamma }\left(t_{i-1},t_{i-1}\right)=0.$$

The use of EKF for FDI problem can be based on the estimation of DS errors. However, in this case, it is necessary to associate each DS state with its own equation of the form (3) and its own filter. There is a need to create a bank of filters [2], which is difficult to implement in practice. Therefore, it is advisable to form FDI algorithms based on the serviceable states of the DS. Then, the diagnostic parameters will reflect the deviation of the actual state of the DS from the serviceable one.

2. Fault detection by the ${\mathrm{\chi }}^2$ Chi-Square criterion

The functioning of the EKF is based on the use of a vector of residuals (see Figure 1)

$${\mathrm{\nu }}_i=Z_i-{\widehat{\mathrm{Z}}}_i={\left[{\mathrm{\nu }}_{1\left(i\right)}{\mathrm{\nu }}_{2\left(i\right)}\dots {\mathrm{\nu }}_{j\left(i\right)}\dots {\mathrm{\nu }}_{l\left(i\right)}\right]}^{\mathrm{T}},$$E7

where ${\widehat{Z}}_i=H_i{\widehat{x}}_{i/i-1}$ and ${\widehat{x}}_{i/i-1}={\mathrm{\Phi }}_i{\widehat{x}}_{i-1/i-1}$ is the predicted values of estimates of the DS errors at the instant $t=t_i$ of time after $Z_{i-1}$ observations have been processed.

The statistical properties of the residuals make it possible to form diagnostic parameters based on them. It is known [2] that the serviceable state of the DS and acceptable error values correspond to a vector of residuals that has a Gaussian distribution with zero mathematical expectation and a covariance matrix ${\mathrm{\alpha }}_i$, that is,

$${\mathrm{\nu }}_i\in N\left(0,{\mathrm{\alpha }}_i^2\right).$$E8

Taking into account the principle of orthogonality of optimal estimates [1] $E\left[e_i{\mathrm{\vartheta }}_i^{\mathrm{T}}\right]=0$, it can be shown that

$${\mathrm{\alpha }}_i^2=E\left[{\mathrm{\nu }}_i{\mathrm{\nu }}_i^T\right]=H_i{P}_{i/i-1}{H}_i^T+R_i,$$E9

where $P_{i/i-1}=E\left[e_{i/i-1}{e}_{i/i-1}^T\right]={\mathrm{\Phi }}_i{P}_{i-1/i-1}{\mathrm{\Phi }}_i^T+{\mathrm{\Gamma }}_i{Q}_{i-1}{\mathrm{\Gamma }}_i^T$ is the predicted covariance matrix; $e_{i/i-1}=x_i-{\widehat{x}}_{i/i-1}$.

In traditional EKF, all elements of the ${\mathrm{\nu }}_i$ vector are processed simultaneously. Taking into account such processing, to monitor the state of the DS, it is necessary to check the l-dimensional vector of residuals for Gaussianity. In practice, the solution of this problem is a matter of some difficulty. Therefore, using the residual vector, more compact diagnostic parameters are formed. One of these parameters is formed by convolving the vector (7) and the covariance matrix (9)

$$J_i={\mathrm{\nu }}_i^T{\mathrm{\alpha }}_i^{-2}{\mathrm{\nu }}_i.$$E10

In the quadratic form (10), the elements of the ${\mathrm{\alpha }}_i$ matrix are the normalizing coefficients that take into account the required statistical characteristics of the residuals. It can be shown [2] that if the ${\mathrm{\nu }}_i$ vector has a Gaussian distribution, then the quadratic form (10) is distributed according to the ${\mathrm{\chi }}^2$ law with l degrees of freedom, namely,

$$J_i\in {\mathrm{\chi }}^2\left(l,2l\right),$$E11

that is, the l dimension of the ${\mathrm{\nu }}_i$ vector is equal to the mathematical expectation of the $J_i$ parameter and to one-half of its variance.

Relation (11) reflects the necessary condition for the absence of violations in the DS. Quadratic form (10) includes all l components of the ${\mathrm{\nu }}_i$ vector, so it can be considered as a DS generalized state parameter. It is believed that if the $J_i$ parameter is within acceptable limits, then there are no violations in the DS. The tolerances for the $J_i$ parameter are formed taking into account the numerical characteristics of the ${\mathrm{\chi }}^2$ distribution and a given level of significance for the goodness-of-fit criterion. Such conditions can be performed taking into account the properties of the a(l) quantile for the ${\mathrm{\chi }}^2$ distribution, namely,

$$P\left\{J_i>J_{T,a\left(l\right)}\right\},$$E12

where $P\left\{\dots \right\}$ is a probability measure; $J_{T,a\left(l\right)}$ is a tabular value of the $J_i$ parameter for the quantile a(l) with the number l of degrees of freedom; a(l) is the quantile of a order, $a\in \left(0,1\right)$.

An analysis of the tabular data shows that quantile (12) with highly accuracy reflects the three mean square deviations rule (the 3σ rule). This rule finds application in the processing of scalar parameters $\mathrm{\xi }$ distributed according to the Gauss law and boils down to the following. The probabilities of deviations of such parameters from the mathematical expectation are characterized by the following values [3]:

$$P\left\{\left|\mathrm{\xi }‐E\left[\mathrm{\xi }\right]\right|\ge k\mathrm{\sigma }\right\}=\left\{\begin{array}{c}0.3173\dots ,k=1;\\ 0.0455\dots ,k=2;\\ 0.0027\dots ,k=3.\end{array}\right.$$E13

For k = 3, relation (13) reflects the 3σ rule for the Gaussian distribution law. Taking into account the probabilistic characteristics (13) and the 0.02(l) quantile of [4], it can be argued that, with a confidence probability of 0.98, the necessary condition for the $J_i$ parameter to belong to the ${\mathrm{\chi }}^2$ distribution is associated with the fulfillment of the following inequality:

$$J_i\le {\mathrm{\gamma }}_l^2=E\left[J_i\right]+3\sqrt{V\left[J_i\right]}=l+3\sqrt{2l},$$E14

where V […] is an operator for variance.

Taking into account the ${\mathrm{\gamma }}_l^2$ tolerance, the monitoring of DS using the diagnostic parameter $J_i$ for the ${\mathrm{\chi }}^2$ criterion is reduced to checking the following conditions:

$$\left.\begin{array}{l}\text{if}{J}_i\le {\mathrm{\gamma }}_l^2,\text{there}\mathrm{are}\mathrm{no}\text{faults in the}\mathrm{DS};\\ \text{if}{J}_i>{\mathrm{\gamma }}_l^2,\text{there}\mathrm{are}\text{faults in the}\mathrm{DS}.\end{array}\right\}$$E15

Monitoring using parameter (10) allows you to determine the state of the DS as a whole without diagnosing which element of the observation vector the violation is associated with. This approach is traditional [5, 6] for DS monitoring. In practice, there is a need to diagnose DS with a depth up to the element of the observation vector.

3. Fault diagnosis by the ${\mathrm{\chi }}^2$ Chi-Square criterion

The problem of DS diagnosis with a depth of up to elements of the vector $Z_{i-1}$ can be solved if the observation errors ${\mathrm{\vartheta }}_i$ are statistically independent (uncorrelated), that is, if the matrix $R_i$ in relation (9) is diagonal. If the observation errors are mutually correlated, their preliminary decomposition can be performed. Such a decomposition can be performed taking into account the structure of the covariance $R_i$. As the specified matrix is square, positive definite and symmetric, it can be represented in the following form [7]:

$$R_i=S_i{D}_i{S}_i^T,$$E16

where $S_i$ is an upper triangular matrix with unit diagonal elements and $D_i$ is a diagonal matrix.

Then, the vector of observations $S_i^{-1}{Z}_i$ will have elements with mutually uncorrelated noise. Indeed, the following model corresponds to such observations:

$$S_i^{-1}{Z}_i=S_i^{-1}{H}_i{x}_i+S_i^{-1}{\mathrm{\vartheta }}_i$$E17

or

$${\tilde{Z}}_i={\tilde{H}}_i{x}_i+{\tilde{\mathrm{\vartheta }}}_i,$$E18

where

$${\tilde{Z}}_i=S_i^{-1}{Z}_i;{\tilde{H}}_i=S_i^{-1}{H}_i;{\tilde{\mathrm{\vartheta }}}_i=S_i^{-1}{\mathrm{\vartheta }}_i.$$

The covariance matrix has the form:

$${\tilde{R}}_i=E\left[{\tilde{\mathrm{\vartheta }}}_i{\tilde{\mathrm{\vartheta }}}_k^T\right]=S_i^{-1}E\left[{\mathrm{\vartheta }}_i{\mathrm{\vartheta }}_k^T\right]S_i^{-T}=S_i^{-1}{R}_i{S}_i^{-T}=S_i^{-1}{S}_i{D}_i{S}_i^T{S}_i^{-T}=D_i,$$E19

where $D_i=\mathit{diag}\left\{D_{1\left(i\right)},\dots D_{j\left(i\right)},\dots D_{m\left(i\right)}\right\}$; $S_i^{-T}={\left(S_i^{-1}\right)}^T$.

Expression (19) shows that ${\tilde{Z}}_i$ observations have uncorrelated errors. Next, we will assume that the matrix $R_i$ is diagonal. Then, such observations can be processed using a sequential modification of the EKF [8].

Prediction:

$$m_0={\widehat{x}}_{i/i-1}={\mathrm{\Phi }}_i{\widehat{x}}_{i-1/i-1};$$E20

$$M_0=P_{i/i-1}={\mathrm{Ф}}_i{P}_{i-1/i-1}{\mathrm{\Phi }}_i^T+{\mathrm{\Gamma }}_i{Q}_{i-1}{\mathrm{\Gamma }}_i^T.$$E21

Updating:

$${\mathrm{\nu }}_j=z_j-{Н}_j{m}_{j-1};$$E22

$${\mathrm{\alpha }}_j^2=H_j{M}_{j-1}{H}_j^T+R_j;$$E23

$$K_j={М}_{j-1}{H}_j^T/{\mathrm{\alpha }}_j^2;$$E24

$$M_j=\left(K_j{H}_j-I\right)M_{j-1}{\left(K_j{H}_j-I\right)}^T+K_j{R}_j{K}_j^T;j=\overline{1,l};$$E25

$${\widehat{x}}_{i/i}=m_l;P_{i/i}=M_l,$$E26

where ${\widehat{x}}_{i/i-1}$ and ${\widehat{x}}_{i/i}$, respectively, are the predicted and updated estimates of the DS error vector at the ith time point, obtained after processing observations $Z_{i-1}$ and $Z_i$;

$P_{i/i-1}$ and $P_{i/i}$ covariance matrices of these estimates; $m_j$ and $M_j$, respectively, estimates of the DS error vector and their covariance matrix obtained after processing the jth observation; $H_i={\left[H_1^T\dots H_j^T\dots H_l^T\right]}^T$ is the matrix for the relation of observed parameters with the DS error vector; $H_j$ is the row vector for the relation of the $z_j$ observation with the DS errors vector; $R_j$ is the variance of the $z_j$ observation error; $Q_i$ is a matrix of weighting coefficients characterizing the noise intensities.

Unlike the traditional EKF, in algorithm (22)–(26), when calculating the amplification coefficient (24), it is not necessary to invert the covariance matrix (9) for the residual vector (7). Such inversion is replaced by l division operations in Eq. (24). In addition, based on the presented algorithm, it seems possible to form statistical criteria of agreement for scalar diagnostic parameters.

When processing the residual vector (7) element-by-element, algorithm (22)–(26) can be supplemented with procedures for analyzing the DS state for each of the l observations. For example, the checking of DS on the jth observation can be performed using a normalized residual ${\mathrm{\beta }}_j={\mathrm{\nu }}_j/{\mathrm{\alpha }}_j$. The statistical properties of such residuals make it possible to form scalar goodness-of-fit criteria. By analogy with the diagnostic parameter (10), for l = 1, the necessary condition for the absence of violations in the DS for each of the observations will have the form:

$${\mathrm{\beta }}_j^2\in {\mathrm{\chi }}^2\left(1,2\right)$$E27

or by the rule of $3\mathrm{\sigma }$, for the quantile a(1) = 0.02,

$${\mathrm{\beta }}_j^2\le {\mathrm{\gamma }}_1^2=E\left[{\mathrm{\beta }}_j^2\right]+3\sqrt{V\left[{\mathrm{\beta }}_j^2\right]}=1+3\sqrt{2}\approx 5.2.$$E28

Based on reference data [4], quantiles $\tilde{a}\left(1\right)=0.01$ and $\overline{a}\left(1\right)=0.001$ can be matched to the extended ${\tilde{\mathrm{\gamma }}}_1^2\approx 7.6$ and ${\overline{\mathrm{\gamma }}}_1^2\approx 11.8$ tolerance values. However, when the tolerances are expanded, the probability of fault detection errors associated with a “false alarm” increases.

Taking into account the ${\mathrm{\gamma }}_1^2$ tolerance, the monitoring of DS using the diagnostic parameter ${\mathrm{\beta }}_j^2$ for the ${\mathrm{\chi }}^2$ criterion is reduced to checking the following conditions:

$$\left.\begin{array}{l}\text{if}{\mathrm{\beta }}_j^2\le {\mathrm{\gamma }}_1^2,\text{there}\mathrm{are}\mathrm{no}\text{faults in the}j\mathrm{th}\text{channel of the}\mathrm{DS};\\ \text{if}{\mathrm{\beta }}_j^2>{\mathrm{\gamma }}_1^2,\text{there}\mathrm{are}\text{faults in the}j\mathrm{th}\text{channel of the}\mathrm{DS}.\end{array}\right\}$$E29

Thus, the decomposition of the residual vector allows the formation of scalar goodness-of-fit criteria. The use of such criteria makes it possible to increase the depth of diagnosis.

4. Fault diagnosis by the ${\mathrm{\vartheta }}^2$ fisher criterion

The use of the ${\mathrm{\chi }}^2$ criterion makes it possible to detect current DS faults in real time. In practice, there is also a need to analyze the DS operation over a certain period of time.

Appropriate diagnostic parameters can be formed based on retrospective data. Technology (22)–(26) for sequential processing of observations allows the formation of scalar diagnostic parameters based on a sample of residuals on a sliding time interval. To do this, you can use the ergodic properties of the estimates obtained by the Kalman filter, the parameters of which are a priori adjusted for the proper operation of the DS. These parameters are the variances of the innovation sequences (22) in each of the l observation channels. The predicted value of the variance of the residual (22) at the ith moment of time is determined by the ratio (23), and its estimate ${\widehat{\mathrm{\alpha }}}_{j\left(i\right)}$ is based on a sample of N values, that is,

$${\widehat{\mathrm{\alpha }}}_{j\left(i\right)}^2=\frac{1}{N-1}\sum _{k=i-N+1}^i{\left[{\mathrm{\nu }}_{j\left(k\right)}-{\overline{\mathrm{\nu }}}_{j\left(i\right)}\right]}^2;j=\overline{1,l},$$E30

where ${\overline{\mathrm{\nu }}}_{j\left(i\right)}=\frac{1}{N}{\sum }_{k=i-N+1}^i{\mathrm{\nu }}_{j\left(k\right)}$ is an estimate of the mathematical expectation of the jth residual at the ith instant of time. Taking into account condition (8), the absence of faults can be associated with the desired value of the mathematical expectation of the residuals, namely, ${\overline{\mathrm{\nu }}}_{j\left(i\right)}=0$. As a parameter that characterizes the DS state on the time interval $T=\left[t_{i‐N+1},t_i\right]$, the ratio of the actual variance ${\widehat{\mathrm{\alpha }}}_j^2$ to its predicted value ${\mathrm{\alpha }}_j^2$ can be accepted, that is,

$$F_j={\widehat{\mathrm{\alpha }}}_j^2/{\mathrm{\alpha }}_j^2.$$E31

It is known [4] that when the condition (32) presented below is satisfied, parameter (31) has ${\mathrm{\vartheta }}^2$ or Fisher distribution, namely,

$$F_j\in {\mathrm{\vartheta }}^2\left(b,c\right),$$E32

$$\text{where}b=N/\left(N–2\right)\text{and}c=4N\left(N-1\right)/\left[{\left(N-2\right)}^2\left(N-4\right)\right]$$E33

are tabular values, respectively, of the mathematical expectation and variance for the $F_j$ parameter.

Taking into account the 3σ rule and a = 0.02 quantile, condition (32) can be represented as

$$F_j\le {\mathrm{\eta }}_1^2=E\left[F_j\right]+3\sqrt{V\left[F_j\right]}=b+3\sqrt{c}.$$E34

When the tolerance ${\mathrm{\eta }}_1^2$ is taken into account, the DS diagnosis according to the ${\mathrm{\vartheta }}^2$ criterion using the $F_j$ parameter is reduced to checking the following conditions:

$$\text{if}{F}_j\le {\mathrm{\eta }}_1^2,$$E35

there were no faults in the DS for the jth channel of observations on the time interval $T=\left[t_{i-N+1},t_i\right]$;

$$\text{if}{F}_j>{\mathrm{\eta }}_1^2,$$E36

there were faults in the DS for the jth channel of observations on the time interval T.

The procedures (35) and (36) complements the test (29) in order to increase the reliability of the diagnosis.

5. Fault diagnosis by the ${\mathrm{\chi }}^2$ and the ${\mathrm{\vartheta }}^2$ combined criterion

When diagnosing DS, the problem of detecting short-term outliers against a background of failures arises. If short-term outliers are recognized, then they can be parried without affecting the operation of the DS. The problem of recognizing the type of violations can be solved by combining the ${\mathrm{\chi }}^2$ and the ${\mathrm{\vartheta }}^2$ goodness-of-fit tests. Indeed, DS diagnosis based on the AA criterion makes it possible to immediately detect all violations, both random short-term outliers and failures. Indeed, DS diagnosis based on the ${\mathrm{\chi }}^2$ criterion makes it possible to immediately detect all violations, both random short-term outliers and failures. At the same time, DS diagnosis based on the ${\mathrm{\vartheta }}^2$ criterion is implemented using a set of residuals on a sliding time interval. Random short-term outliers that may appear on a sliding time interval when calculating the $F_j$ parameter are averaged and, as a rule, do not affect to the diagnosis results of the ${\mathrm{\vartheta }}^2$ test. On the other hand, the appearance of gradual and sudden failures, which are characterized by constant deviations of residuals from their nominal values, lead to the $F_j$ parameter going beyond the tolerance limits. Thus, if a violation in the jth observation was detected according to both criteria, then most likely a permanent failure occurred in the DS, and if only according to the ${\mathrm{\chi }}^2$ criterion, then it was a random short-term outlier.

From the considered properties of the combined criterion of consent, the following algorithm for detecting and parrying violations in the observed DS is as follows:

• If there is no discrepancy between the diagnostic parameters relative to the tolerances, the residuals are processed using EKF.

• Anomalous observation signals detected by the ${\mathrm{\chi }}^2$ criterion are excluded from processing or processed with certain robust confidence coefficients.

• Violations detected by both criteria ${\mathrm{\chi }}^2$ and ${\mathrm{\vartheta }}^2$ are countered by connecting a redundant channel.

The procedures presented above allow for the diagnosis of DS with a depth of up to the element of the observations vector. At the same time, in practice, it becomes necessary to perform diagnosis with a depth up to the element of the DS errors vector.

6. Fault diagnosis based on joint filtering and smoothing procedures

The diagnosis of DS with a depth up to the elements of the errors vector can be implemented on the basis of joint procedures for filtering observations in “forward” time and smoothing estimates in “backward” time. By analogy with expression (10), the diagnostic parameter, which reflects the deviations of error vector estimates from their predicted root-mean-square values, can be represented by a quadratic form:

$$J_i={\mathrm{\nu }}_{i/N}^T\mathrm{\Delta }{P}_i^{-1}{\mathrm{\nu }}_{i/N},$$E37

where ${\mathrm{\nu }}_{i/N}={\mathrm{\delta }}_{f\left(i\right)}-{\mathrm{\delta }}_{\mathrm{s}\left(i\right)}={\mathrm{\Phi }}_{i+1}^{-1}{\widehat{x}}_{i+1/N}-{\widehat{x}}_{i/i}$;

${\mathrm{\delta }}_{f\left(i\right)}=x_i-{\widehat{x}}_{i/i}$; ${\mathrm{\delta }}_{s\left(i\right)}={\mathrm{\Phi }}_{i+1}^{-1}\left(x_{i+1}-{\widehat{x}}_{i+1/N}\right)$;

$\mathrm{\Delta }{P}_i=P_{i/i}+{\mathrm{\Phi }}_i^{-1}{P}_{i+1/N}{{\mathrm{\Phi }}_i}^{-T};$

${\widehat{x}}_{i/i}$ and ${\widehat{x}}_{i/N}$ are estimates of the DS error vector at the ith moment of time, obtained, respectively, at the filtering and smoothing stage; $P_{i/i},P_{i/N}$ are the covariance matrices of the above estimates; ${\mathrm{Ф}}^{-Т}={\left({\mathrm{Ф}}^{‐1}\right)}^{Т}$. ${\mathrm{\Phi }}_i^{-1}$ is the inverse transition matrix, which is determined from the solution of the differential equation:

$$d{\mathrm{\Phi }}^{-1}\left(t,t_{i-1}\right)/\mathit{dt}=-{\mathrm{Ф}}^{-1}\left(t,t_{i-1}\right)A\left(t\right);{\mathrm{\Phi }}^{-1}\left(t_{i-1},t_{i-1}\right)=I.$$E38

Stable smoothing (∣δ∣ < 3σ), which reflects the absence of anomalous observations signals, is characterized by the following distribution of the residual ${\mathrm{\nu }}_{i/N}$ and the quadratic form $J_{s\left(i\right)}$:

$${\mathrm{\nu }}_{i/N}\in N\left(0;\mathrm{\Delta }{P}_i\right);J_{s\left(i\right)}\in {\mathrm{\chi }}^2\left(n;2n\right),$$E39

where n is the dimension of the DS error vector.

Taking into account the properties of the ${\mathrm{\chi }}^2$ distribution and the 3σ rule, it is possible to formulate the necessary conditions for the absence of violations for the DS in general and for the jth element of the error vector in particular

$$J_{s\left(i\right)}\le n+3\sqrt{2n};$$E40

$$J_{s\left(i/j\right)}=J_{s\left(i/j-1\right)}+{\tilde{\mathrm{\nu }}}_{i/N\left(j\right)}^2/\mathrm{\Delta }{D}_{i\left(j\right)}\le {\mathrm{\gamma }}_j^2=j+3\sqrt{2j},$$E41

where ${\tilde{\mathrm{\nu }}}_{i/N}=\mathrm{\Delta }{U}_i^{-1}{\mathrm{\nu }}_{i/N}$; $j=\overline{1,n}$; $\mathrm{\Delta }{U}_i^{-1};\mathrm{\Delta }{D}_i^{-1}$ are an upper triangular matrix with identity diagonal and a diagonal matrix, respectively. These matrices are obtained by the following orthogonal transformation:

$$\mathrm{\Delta }{P}_i^{-1}=\mathrm{\Delta }{U}_i^{-T}\mathrm{\Delta }{D}_i^{-1}\mathrm{\Delta }{U}_i^{-1};$$E42

$\mathrm{\Delta }{D}_{i\left(j\right)}^{-1}$ is the jth element of the diagonal matrix $\mathrm{\Delta }{D}_i^{-1}$.

Taking into account decomposition (42) and the statistical properties of the Fisher distribution [4]:

$$F_{s\left(j\right)}={\widehat{\mathrm{\alpha }}}_{\left(i/j\right)}/\mathrm{\Delta }{D}_{\left(i/j\right)}\in {\mathrm{\vartheta }}^2\left(b,c\right),$$E43

it is possible to formulate the necessary condition for the operable DS state (“absence” of failures) for the jth element of the error vector, that is,

$$F_{s\left(j\right)}\le {\mathrm{\eta }}_j^2=b+3\sqrt{c},$$E44

where ${\widehat{\mathrm{\alpha }}}_{\left(i/j\right)}$ is an estimate of the variance of the residual ${\tilde{\mathrm{\nu }}}_{i/N\left(j\right)}$ on a sliding time interval; b and c are tabulated values (33) of the mathematical expectation and variance for the parameter $F_j$.

Diagnosis based on joint procedures for optimal filtering and smoothing assumes the relationship between prediction and interpolation models of the DS error vector:

$$x_{i+1}={\mathrm{\Phi }}_{i+1}{x}_i+{\mathrm{\Gamma }}_{i+1}{\mathrm{\xi }}_i\text{prediction one step ahead at}i=\overline{i_0,i_f};$$E45

$$x_i={{\mathrm{\Phi }}^{-1}}_{i+1}\left(x_{i+1}-{\mathrm{\Gamma }}_{i+1}{\mathrm{\xi }}_i\right)\text{interpolation one step back at}i=\overline{i_f,i_0}.$$E46

It should be noted, however, that the traditional RTS (Rauch-Tung-Striebel) algorithm [9] smooths out the predicted estimates ${\widehat{x}}_{i/i-1}$, the reliability of which is significantly lower than the updated estimates ${\widehat{x}}_{i/i}$.

The value ${\widehat{x}}_{\left(i/i\right)N}$ for the smoothed estimate of the $x_i$ error vector can be found by solving the following optimization problem using the least squares method

$${\widehat{x}}_{\left(i/i\right)N}=\underset{x_i}{argmin}{J}_{i/s},$$E47

$$\text{where}\begin{array}{l}{J}_{i/s}=0.5\{\underset{a}{\underbrace{{\left[x_i-{\widehat{x}}_{\left(i/i+1\right)N}\right]}^T{P}_{\left(i/i+1\right)N}^{-1}\left[x_i-{\widehat{x}}_{\left(i/i+1\right)N}\right]}}+\\ +\underset{b}{\underbrace{{\left[x_i-{\widehat{x}}_{i/i}\right]}^T{P}_{i/i}^{-1}\left[x_i-{\widehat{x}}_{i/i}\right]\}}}\end{array};$$E48

${\widehat{x}}_{\left(i/i+1\right)N}$; $P_{\left(i/i+1\right)N}$ is an estimate of the error vector interpolated at the ith moment of time and its covariance matrix.

The quadratic form (48) combines components “a” and “b”, which characterize interpolation and filtering errors at the ith moment of time. It should be noted that component “b” can also be formed for predicted estimates. However, in this case, smoothing errors will significantly depend on the frequency of registration these estimates.

Applying the rule of differentiation with respect to the parameter for the quadratic form (48), we can obtain the following solution to problem (47)

$${\left.\frac{\partial J_i}{\partial x_i^T}\right|}_{x_i={\widehat{x}}_{\left(i/i\right)N}}=P_{\left(i/i+1\right)N}^{-1}\left[{\widehat{x}}_{\left(i/i\right)N}-{\widehat{x}}_{\left(i/i+1\right)N}\right]+P_{i/i}^{-1}\left[{\widehat{x}}_{\left(i/i\right)N}-{\widehat{x}}_{i/i}\right]=0.$$E49

From (Eq. (49)), an expression can be obtained for a smoothed estimate of the error vector

$${\widehat{x}}_{\left(i/i\right)N}={\widehat{x}}_{\left(i/i+1\right)N}-P_{\left(i/i\right)N}{P}_{i/i}^{-1}\left[{\widehat{x}}_{\left(i/i+1\right)N}-{\widehat{x}}_{i/i}\right],$$E50

$$\text{where}{P}_{\left(i/i\right)N}={{\left[P_{\left(i/i+1\right)N}^{-1}+P_{i/i}^{-1}\right]}^{-1}}_{;}$$E51

$$K_{i/N}=P_{\left(i/i\right)N}{P}_{i/i}^{-1}\text{is the optimal amplification coefficient}.$$E52

Applying the matrix inversion lemma [10] to expression (51), it can be represented in expanded form:

$$P_{\left(i/i\right)N}=P_{\left(i/i+1\right)N}-P_{\left(i/i+1\right)N}{\left[P_{\left(i/i+1\right)N}+P_{i/i}\right]}^{-1}{P}_{\left(i/i+1\right)N}.$$E53

Substituting expression (53) into relation (52), we obtain the following procedure for calculating the optimal amplification coefficient

$$K_{i/N}=P_{\left(i/i+1\right)N}{\left[P_{\left(i/i+1\right)N}+P_{i/i}\right]}^{-1}.$$E54

The expression for the interpolated estimate follows from the ratio (46)

$${\widehat{x}}_{\left(i/i+1\right)N}={\mathrm{Ф}}_i^{-1}{\widehat{x}}_{\left(i+1/i+1\right)N}.$$E55

The difference between relations (46) and (55) determines the equation for interpolation errors

$$e_{\left(i/i+1\right)N}={\mathrm{Ф}}_i^{-1}\left[e_{\left(i+1/i+1\right)N}-{\mathrm{Г}}_i{\mathrm{\xi }}_i\right],$$E56

where $e_{\left(i/i+1\right)N}=x_i-{\widehat{x}}_{\left(i/i+1\right)N}$; $e_{\left(i+1/i+1\right)N}=x_{i+1}-{\widehat{x}}_{\left(i+1/i+1\right)N}$.

The following covariance matrix corresponds to the interpolation error vector

$$P_{\left(i/i+1\right)N}=E\left[e_{\left(i/i+1\right)N}{e}_{\left(i/i+1\right)N}^T\right]={\mathrm{Ф}}_i^{-1}\left[P_{\left(i+1/i+1\right)N}+{\mathrm{Г}}_i{Q}_i{\mathrm{Г}}_i^T\right]{\mathrm{Ф}}_i^{-T}.$$E57

Combining (Eqs. (50), (53)–(55), (57)) into a single structure, we obtain an algorithm for optimal smoothing of error vector estimates

Interpolation:

$${\widehat{x}}_{\left(i/i+1\right)N}={\mathrm{Ф}}_{i+\mathit{1}}^{-1}{\widehat{x}}_{\left(i+1/i+1\right)N};$$E58

$${\tilde{P}}_{\left(i/i+1\right)N}=P_{\left(i+1/i+1\right)N}+{\mathrm{Г}}_{i+1}{Q}_i{\mathrm{Г}}_{i+1}^T.$$E59

Updating:

$${\mathrm{\nu }}_{i/N}={\widehat{x}}_{\left(i/i+1\right)N}-{\widehat{x}}_{i/i};$$E60

$$P_{\left(i/i\right)N}^{-1}={\mathrm{Ф}}_{i+1}^T{\tilde{P}}_{\left(i/i+1\right)N}^{-1}{\mathrm{Ф}}_{i+1}+P_{i/i}^{-1};$$E61

$$K_{i/N}=P_{\left(i/i\right)N}{P}_{i/i}^{-1};$$E62

$${\widehat{x}}_{\left(i/i\right)N}={\widehat{x}}_{\left(i/i+1\right)N}-K_{i/N}{\mathrm{\nu }}_{i/N},$$E63

To implement the algorithm (58)–(63), it is necessary to register estimates ${\widehat{x}}_{i/i}$of the DS error vector and their covariance matrices $P_{i/i}^{-1}$ obtained by the EKF in the process of functional monitoring of the DS. The matrix $P_{i/i}^{-1}$ is determined from the solution in “forward” time according to the following equations:

$$P_{i/i}^{-1}=P_{i/i-1}^{-1}+H_i^T{R}_i^{-1}{H}_i;$$E64

$$P_{i/i-1}^{-1}={\left({\mathrm{\Phi }}_i{P}_{i-1/i-1}{\mathrm{\Phi }}_i^T+{\mathrm{\Gamma }}_i{Q}_{i-1}{\mathrm{\Gamma }}_i^T\right)}^{-1}.$$E65

The implementation of (Eqs. (58)–(65)) is associated with computational difficulties due to the need to determine the inverse covariance matrices $P_{i/i}^{-1}$ and $P_{\left(i/i\right)N}^{-1}$. These difficulties can be significantly reduced with the sequential implementation of smoothing algorithms. With this implementation, the inversion of covariance matrices of n×n dimension is replaced by a sequence of division operations. Such a replacement is based on the transformation of (Eq. (61), (64), (65)) according to the matrix inversion lemma [10]. For (Eq. (65)) it will have the form:

$$P_{i/i-1}^{-1}={\tilde{P}}_{i/i-1}^{-1}-{\tilde{P}}_{i/i-1}^{-1}{\mathrm{\Gamma }}_i{\left({\mathrm{\Gamma }}_i^T{\tilde{P}}_{i/i-1}^{-1}{\mathrm{\Gamma }}_i+Q_{i-1/i-1}^{-1}\right)}^{-1}{\mathrm{\Gamma }}_i^T{\tilde{P}}_{i/i-1}^{-1}={\left({\tilde{P}}_{i/i-1}+\sum _{j=i}^r{\mathrm{\Gamma }}_j{\mathrm{\Gamma }}_j^T{Q}_{\mathit{jj}}\right)}^{-1},$$E66

where ${\mathrm{\Gamma }}_j$ is the jth column of matrix ${\mathrm{\Gamma }}_i$, which has n×r dimension; $Q_{\mathit{jj}}$ is the jth element of the diagonal matrix $Q_i$, which has r×r dimension.

Hence, the following procedure for the sequential implementation of Eq. (66) is valid

$$M_0^{-1}={\mathrm{\Phi }}_i^{-T}{P}_{i-1/i-1}^{-1}{\mathrm{\Phi }}_i^{-1};$$E67

$$K_j=M_{j-1}^{-1}{\mathrm{\Gamma }}_j/\left({\mathrm{\Gamma }}_j^T{M}_{j-1}^{-1}{\mathrm{\Gamma }}_j+Q_{\mathit{jj}}^{-1}\right);$$E68

$$M_j^{-1}=\left(K_j{\mathrm{\Gamma }}_j^T-I\right)M_{j-1}^{-1}{\left(K_j{\mathrm{\Gamma }}_j^T-I\right)}^T+K_j{Q}_{\mathit{jj}}^{-1}{K}_j^T;$$E69

$$P_{i/i-1}^{-1}=M_r^{-1};j=\overline{1,r}.$$E70

Transformations similar to (66)–(70) will also be valid for Eq. (61), if the covariance matrix of filtering errors is presented in the following form

$$P_{i/i}^{-1}=U_{i/i}^{-T}{D}_{i/i}^{-1}{U}_{i/i}^{-1},$$E71

where $U_{i/i}$ is an upper triangular matrix with unit diagonal elements; $D_{i/i}$ is a diagonal matrix. These matrices form the basis of the U-D filter [11].

Taking into account decomposition (71) and preliminary decomposition of residuals ${\tilde{\mathrm{\nu }}}_{i/N}=U_i^{-1}{\mathrm{\nu }}_{i/N}$, a sequential modification of the smoothing and diagnostic algorithm can be formed.

Interpolation:

$$m_0={\tilde{x}}_{\left(i/i+1\right)N}={\mathrm{\Phi }}_{i+1}^{-1}{\widehat{x}}_{\left(i+1/i+1\right)N};$$E72

$${\mathrm{\nu }}_{i/N}={\tilde{x}}_{\left(i/i+1\right)N}-{\widehat{x}}_{i/i};$$E73

$$\mathit{\text{MWGS}}\left\{\begin{array}{c}{\overline{W}}_0=\left[{\mathrm{\Phi }}_{i+1}^{-1}{U}_{i+1/N}⋮{\mathrm{\Phi }}_{i+1}^{-1}{\mathrm{\Gamma }}_{i+1}\right]\\ {\overline{D}}_0=\mathit{diag}\left(D_{i+1/N},Q_i\right)\end{array}\right\}\to \begin{array}{c}{\tilde{U}}_0\\ {\tilde{D}}_0\end{array};$$E74

Smoothing:

$$f_j={\overline{U}}_j^{-1}{\tilde{U}}_{j-1};V_j={\tilde{D}}_j{f}_j^T;$$E75

$${\mathrm{\alpha }}_j^2=f_j{V}_j+D_{\mathit{jj}};$$E76

$${\tilde{K}}_j={\tilde{U}}_{j-1}{V}_j/{\mathrm{\alpha }}_j^2;$$E77

$$\mathit{\text{MWGS}}\left\{\begin{array}{c}{\overline{W}}_j=\left[\left({\tilde{K}}_j{f}_j-{\tilde{U}}_{j-1}\right)⋮{\tilde{K}}_j\right]\\ {\overline{D}}_j=\mathit{diag}\left({\tilde{D}}_{j-1};D_{\mathit{jj}}\right)\end{array}\right\}\to \begin{array}{c}{\tilde{U}}_j\\ {\tilde{D}}_j\end{array};$$E78

$$K_j={\tilde{U}}_j{\tilde{D}}_j{\tilde{U}}_j^T{\overline{U}}_j^{-T}{D}_j^{-1};$$E79

$${\tilde{\mathrm{\nu }}}_{i/N\left(j\right)}={\overline{U}}_j^{-1}{\mathrm{\nu }}_{i/N};$$E80

$$m_j=m_{j-1}-K_j{\tilde{\nu }}_j;j=\overline{1,n};$$E81

$${\widehat{x}}_{i/N}=m_n;U_{\left(i/i\right)N}={\tilde{U}}_n;D_{\left(i/i\right)N}={\tilde{D}}_n,$$E82

where ${\overline{U}}_j^{-1}$ is the jth row of the matrix $U_{i/i}^{-1}$; $D_j^{-1}$ is the jth element of the diagonal matrix $D_{i/i}^{-1}$; MWGS is the modified weighted Gram-Schmidt procedure [11], intended to transform the aggregate of matrices ${\overline{W}}_j$ and ${\overline{D}}_j$, which are $n\times \left(n+r\right)$ matrix and $\left(n+r\right)\times \left(n+r\right)$ matrix, respectively, into the aggregate of the $n\times n$ an upper triangular matrix ${\tilde{U}}_j$ with identity diagonal and a diagonal matrix ${\tilde{D}}_j$. The algorithm for implementing the MWGS procedure has the form

$$d_m=\sum _{j=1}^{n+r}{\overline{W}}_{\mathit{mj}}{\overline{D}}_j{\overline{W}}_{\mathit{mj}},m=\overline{n,2};$$E83

$$\left.\begin{array}{l}{u}_{\mathit{km}}=\sum _{j=1}^{n+r}{\overline{W}}_{\mathit{mj}}{\overline{D}}_j{\overline{W}}_{\mathit{mj}}/d_m,\\ {\overline{W}}_{\mathit{kj}}≔{\overline{W}}_{\mathit{kj}}-u_{\mathit{km}}{\overline{W}}_{\mathit{mj}},j=\overline{1,n+r}\end{array}\right\}k=\overline{1,m-1};$$E84

$$d_1=\sum _{j=1}^{n+r}{\overline{W}}_{1j}{\overline{D}}_j{\overline{W}}_{1j};$$E85

$u_{\mathit{km}}$ and $d_j$ are the elements of the required matrix U and D.

To implement the U-D smoothing algorithm, it is necessary to determine the $U_{i/i}^{-1}$, $D_{i/i}^{-1}$ components of the a posteriori covariance matrix of estimation errors $P_{i/i}^{-1}$. The traditional U-D algorithm [11] for observation filtering, which includes additional procedures necessary to smoothing estimates of the DS error vector, is presented below.

Prediction for filtering:

$$m_0=x_{i/i-1}={\mathrm{\Phi }}_i{\widehat{x}}_{i-1/i-1};$$E86

$$\mathit{\text{MWGS}}\left\{\begin{array}{c}{\overline{W}}_i=\left[{\mathrm{\Phi }}_i{U}_{i-1/i-1}⋮{\mathrm{\Gamma }}_i\right]\\ {\overline{D}}_i=\mathit{diag}\left(D_{i-1/i-1},Q_{i-1}\right)\end{array}\right\}\to U_0;D_0.$$E87

Prediction for smoothing:

$$U_0^{-T}={\mathrm{\Phi }}_i^{-T}{U}_{i-1/i-1}^{-T};D_0^{-1}≔D_{i-1/i-1}^{-1};$$E88

$$f_j={\mathrm{\Gamma }}_j{U}_{j-1}^{-T};V_j=D_{j-1}^{-1}{f}_j^T;$$E89

$$K_j=U_{j-1}^{-T}{V}_j/\left(f_j{V}_j+Q_{\mathit{jj}}^{-1}\right);$$E90

$$\mathit{\text{MWGSL}}\left\{\begin{array}{c}{\tilde{W}}_j=\left[K_j{f}_j-U_j^{-T}⋮K_j\right]\\ {\tilde{D}}_j=\mathit{diag}\left(D_j^{-1},Q_{\mathit{jj}}^{-1}\right)\end{array}\right\}\to U_j^{-T};D_j^{-1};j=\overline{1,r};$$E91

$$U_0^{-T}≔U_{i/i-1}^{-T}=U_r^{-T};D_0^{-1}≔D_{i/i-1}^{-1}=D_r^{-1}.$$E92

Filtering:

$$f_j=H_j{U}_{j-1};V_j=D_{j-1}{f}_j^T;$$E93

$${\mathrm{\alpha }}_j^2=f_j{V}_j+R_j;$$E94

$$K_j=U_{j-1}{V}_j/{\mathrm{\alpha }}_j^2;$$E95

$$m_j=m_{j-1}+K_j{\mathrm{\nu }}_j;$$E96

$$\mathit{\text{MWGS}}\left\{\begin{array}{c}{\overline{W}}_j=\left[K_j{f}_j-U_{j-1}⋮K_j\right]\\ {\overline{D}}_j=\mathit{diag}\left(D_{j-1},R_j\right)\end{array}\right\}\to U_j;D_j.j=\overline{1,l}.$$E97

Updating for smoothing:

$$\mathit{\text{MWGSL}}\left\{\begin{array}{c}{\tilde{W}}_j=\left[U_{j-1}^{-T}⋮H_j^T\right]\\ {\tilde{D}}_j=\mathit{diag}\left(D_{j-1}^{-1},R_j^{-1}\right)\end{array}\right\}\to U_j^{-T};D_j^{-1};j=\overline{1,l};$$E98

$$U_{i/i}^{-T}=U_l^{-T};D_{i/i}^{-1}=D_l^{-1},$$E99

where MWGSL is a similar procedure to MWGS for forming matrices: the lower triangular with unit diagonal elements $U_j^{-T}$ and diagonal $D_j^{-1}$. Eq. (88)–(92), (98), (99) are redundant with respect to the basic U-D filtering algorithm [11].

Diagnosis in “backward time” is performed when violations detected in “forward time”.

7. Analysis of the results of mathematical modeling of diagnosis procedures based on filtering and smoothing of observations

The Schuler pendulum [12], which includes a gyro-accelerometer (G-A) system moving in a gravitational field along a sphere of R radius, is considered as an object of mathematical modeling. The invariance of the modeled vertical (radius vector) to the movement of the base of the G-A system along the sphere when changing the velocity V and angular position φ is ensured by precession of the gyros. For this purpose, a $\dot{\mathrm{\phi }}=V/R$ signal that is proportional to angular velocity of G-A system motion relative to the sphere is fed to the gyros moment sensor. When the vector $x\left(t\right)={\left[\mathrm{\Delta }V\mathrm{\delta }\mathrm{\Delta }a\mathrm{\Delta \omega }\right]}^{\mathrm{T}}$ of errors of such a system is observed from the velocity signals $Z_V\left(t\right)=V_{\mathrm{G}‐\mathrm{A}}\left(t\right)-V_{\mathrm{SEI}}\left(t\right)$ the parameters of Eqs. (3) and (5) will have the following form

$$A\left(t\right)=\left[\begin{array}{cccc}0& -g& 1& 0\\ 1/R& 0& 0& 1\\ 0& 0& -1/{\mathrm{\tau }}_a& 0\\ 0& 0& 0& -1/{\mathrm{\tau }}_{\mathrm{\omega }}\end{array}\right];$$E100

$$G\left(t\right)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& {\mathrm{\sigma }}_a\sqrt{2/{\mathrm{\tau }}_a}& 0\\ 0& 0& 0& {\mathrm{\sigma }}_{\mathrm{\omega }}\sqrt{2/{\mathrm{\tau }}_{\mathrm{\omega }}}\end{array}\right];$$E101

$$H\left(t\right)={\left[1000\right]}^{\mathrm{T}},$$E102

where g is the acceleration of gravity; δ is an angular error of determining the vertical; $\mathrm{\Delta V}$ is an error of velocity reckoning; $\mathrm{\Delta }a$ is an accelerometer error; $\mathrm{\Delta \omega }$ is the gyro drift; ${\mathrm{\tau }}_a,{\mathrm{\tau }}_{\omega }$ are the correlation time for the accelerometer error and gyro drift, respectively; ${\mathrm{\sigma }}_a,{\mathrm{\sigma }}_{\mathrm{\omega }}$ are mean-square values of errors of the accelerometer and gyro, respectively; $\mathrm{\Delta }\left(\dots \right)$ stands for an error.

In Figures 2 and 3, characteristic results of the studies of an algorithm for G-A system diagnosis from the recorded data are presented. At the 500th second, the accelerometer failure corresponding to the imitation bias was simulated. Such a failure indirectly makes itself evident during filtering for the velocity observation channel when the diagnostic parameter ${\mathrm{\beta }}_V^2$ is above upper tolerance.

Figure 2.

Estimates of the accelerometer error.

Figure 3.

Estimates of the gyro drift.

In post processing the recorded estimates and during the diagnosis by rules (41) and (44), it is possible to determine which of the sensors – an accelerometer or a gyro – has most probably caused a fault.

Figures 2 and 3 show the dynamics of variation of the estimates of both the bias of the $a_x$ accelerometer output signal and the ${\mathrm{\omega }}_x$ gyro drift when velocity observations are processed in “forward time” and when the estimates mentioned are refined in “backward time”. During the diagnosis from the recorded data, a failed accelerometer is isolated if the generalized parameters $J_{{\mathit{Sa}}_j}$ (the ${\mathrm{\chi }}^2$ test) and $F_{{\mathit{Sa}}_j}$ (the ${\mathrm{\vartheta }}^2$ test) are above upper tolerances (see Figure 2). It can also be seen (see Figure 3) that the accelerometer failure has had an insignificant influence on the variation of the generalized parameters $J_{S{\mathrm{\omega }}_j}$ and $F_{S{\mathrm{\omega }}_j}$, which characterize the status of the ${\mathrm{\omega }}_x$ gyro.

Thus, combined processing of observations in “forward” and “backward” time enables us to solve diagnosis problems with a depth of the component of the state vector of a dynamical system.

8. Analysis and parrying of violations by robust filtering

The analysis and parrying of anomalous observations can be performed based on a robust modification of the U-D filter. This modification relies on the application of the influence function ψ, which defines the level of confidence in incoming observations. Such a function can be formed for the normalized residual ${\beta }_j={\nu }_j/{\alpha }_j$, where ${\mathrm{\alpha }}_j$ is a scaling parameter which may be the mean square deviation $\sqrt{R_j}$. For the residual ${\mathrm{\beta }}_j$, one can not only form the robust-likelihood function ρ(β) but also perform the optimization of the estimates

$${\widehat{x}}_i=\underset{x_i}{argmin}\sum _{i=i_0+1}^{i_f}\mathrm{\rho }\left({\mathrm{\beta }}_i\right),$$E103

where $\mathrm{\rho }\left(\mathrm{\beta }\right)=-lnf\left(\mathrm{\beta }\right)$; $f\left(\mathrm{\beta }\right)$ is a probability density function (PDF) [3].

The solution of the problem (116), in view of the constraint $x_i-{\mathrm{Ф}}_i{x}_{i-1}-{\mathrm{Г}}_i{\mathrm{\xi }}_{i-1}=0$, is an algorithm for robust estimation.

Prediction: Eq. (97), (98)

Tuning:

$${\mathrm{\beta }}_j={\mathrm{\nu }}_j/{\mathrm{\alpha }}_j;{\mathrm{\psi }}_j=\mathrm{\psi }\left({\mathrm{\beta }}_j\right);{\mathrm{\psi }}_j^{\prime }={\mathrm{\psi }}^{\prime }\left({\mathrm{\beta }}_j\right);$$E104

Updating:

$$f_j=H_j{U}_{j-1};V_j=D_{j-1}{f}_j^T;$$E105

$${\tilde{\mathrm{\alpha }}}_j^2=f_j{V}_j{\mathrm{\psi }}^{\prime }+{\mathrm{\alpha }}_j^2;$$E106

$$K_j=U_{j-1}{V}_j/{\tilde{\mathrm{\alpha }}}_j^2;$$E107

$$m_j=m_{j-1}+K_j{\alpha }_j{\psi }_j;$$E108

$$\mathit{\text{MWGS}}\left\{\begin{array}{c}{\overline{W}}_j=\left[K_j{f}_j{\mathrm{\psi }}_j^{\prime }-U_{j-1}⋮K_j\right]\\ {\overline{D}}_j=\mathit{diag}\left(D_{j-1},{\mathrm{\alpha }}_j^2{\mathrm{\psi }}_j^{\prime }\right)\end{array}\right\}\to U_j;D_j;$$E109

$$U_{i/i}=U_l;D_{i/i}=D_l;{\widehat{x}}_{i/i}=m_l;j=\overline{1,l},$$E110

where ${\mathrm{\psi }}_j=\mathrm{\psi }\left({\mathrm{\beta }}_j\right)=\left.\frac{\mathrm{\partial \rho }\left(\mathrm{\beta }\right)}{\mathrm{\partial \beta }}\right|\begin{array}{c}\\ \mathrm{\beta }={\mathrm{\beta }}_j\end{array}$ and ${\mathrm{\psi }}_j^{\prime }={\mathrm{\psi }}^{\prime }\left({\mathrm{\beta }}_j\right)=\left.\frac{{\partial }^2\mathrm{\rho }\left(\mathrm{\beta }\right)}{\partial {\mathrm{\beta }}^2}\right|\begin{array}{c}\\ \mathrm{\beta }={\mathrm{\beta }}_j\end{array}$.

The influence function [13, 14] ${\mathrm{\psi }}_j$ and its derivative ${\mathrm{\psi }}_j^{\prime }$ can be formed taking into account “a priori” assumptions about the laws of distribution of the useful signal and noise. The choice of the values of above functions is based on the necessary conditions [15] for the filter to be divergence-free, namely, the diagnostic parameter must have a Gaussian distribution, that is, $\mathrm{\beta }\in N\left(0,1\right)$; in addition, the rule of $3\mathrm{\sigma }$ is fulfilled [3] for the probability P that a random variable having the Gaussian distribution will be on the interval $\left[-3\mathrm{\sigma },3\mathrm{\sigma }\right]$, that is,

$$P\left\{\left|\mathrm{\beta }-E\left[\mathrm{\beta }\right]\right|\ge 3\mathrm{\sigma }\right\}=0.0027.$$E111

For the random variable $\mathrm{\beta }$, the rule (111) can be written as $P\left\{\left|\mathrm{\beta }\right|<3\right\}=0.9973$.

Thus, the correct functioning of the DS can be matched by the condition $\left|\mathrm{\beta }\right|<3$ and also the following values of the functions:

$$f_g\left(\mathrm{\beta }\right)={\left(2\mathrm{\pi }\right)}^{-\frac{1}{2}}exp\left(-0.5{\mathrm{\beta }}^2\right);{\mathrm{\rho }}_g\left(\beta \right)=0.5ln\left(2\mathrm{\pi }\right)+0.5{\mathrm{\beta }}^2;{\mathrm{\psi }}_g\left({\mathrm{\beta }}_j\right)={\mathrm{\beta }}_j;{\mathrm{\psi }}_g^{\prime }\left({\mathrm{\beta }}_j\right)=1.$$E112

The following functions can be made to correspond to such a distribution and to off-design conditions of the filter operation.

$$f_l\left(\mathrm{\beta }\right)=0.5exp\left(-|\mathrm{\beta }|\right);{\mathrm{\rho }}_l\left(\mathrm{\beta }\right)=ln2+\left|\mathrm{\beta }\right|;{\mathrm{\psi }}_l\left({\mathrm{\beta }}_j\right)=\left|1\right|;{\mathrm{\psi }}_l^{\prime }\left({\mathrm{\beta }}_j\right)=0.$$E113

The vagueness of boundaries between anomalous and conditioned signals can be taken into account based on the convolution of the Gauss and Laplace PDF. Such a convolution can be performed using the following moment generating functions (MGFs) [16]:

$$M_{\mathit{lg}}\left(T\right)=M_l\left(T\right)M_g\left(T\right)={\left(1-T^2\right)}^{-1}exp\left(T^2/2\right),$$E114

where M(T) is an MGF; $T$ is generally a complex number; and $K\left(T\right)=lnM\left(T\right)$ is a cumulant generating function.

Relying on the results of [14], the following relations can be shown to hold for the normalized residual ${\mathrm{\beta }}_j$:

$$\mathrm{\psi }\left({\mathrm{\beta }}_j\right)={Т}_0+{\left.K^{\left(3\right)}\left(T\right)/2\right|}_{T=T_0};$$E115

$${\mathrm{\psi }}^{\prime }\left({\mathrm{\beta }}_j\right)=\left[1+\partial K^{\left(3\right)}\left(T\right)/\partial T\right]{\left.\frac{\partial T}{\mathrm{\partial \beta }}\right|}_{T=T_0},$$E116

where ${Т}_0$ is the value of the argument at a saddle point for which the following equality is valid:

$${К}_{Т}^{\prime }\left(Т\right)-{\mathrm{\beta }}_j=0.$$E117

In view of the approximation $ln\left(1-T^2\right)\approx -T^2$ and relations (112)–(117), the parameters ${\mathrm{\psi }}_j$and ${\mathrm{\psi }}_j^{\prime }$ take the form ${\mathrm{\psi }}_{\mathit{lg}}\left({\mathrm{\beta }}_j\right)=\mathrm{\beta }/3$; ${\mathrm{\psi }}_{\mathit{lg}}^{\prime }\left({\mathrm{\beta }}_j\right)=1/3$.

The influence function which reflect the considered assumptions is shown in Figure 4.

Figure 4.

Diagram for the control of an estimating filter with an influence function.

In this figure, a typical influence function is denoted by a dotted line, too. However, the interrelation between such a function and the parameter ${\mathrm{\psi }}_j^{\prime }$ is not obvious.

The need to parry anomalous observations can be illustrated by the results of mathematical modeling of a G-A system, which is described by Eq. (100)–(102). Such results are presented in Figures 5 and 6, where the true $\mathrm{\Delta \delta }={х}_{\mathrm{\delta }}-{\widehat{х}}_{\mathrm{\delta }}$ and root mean square (RMS) error ${\mathrm{\sigma }}_{\mathrm{\delta }}=\sqrt{P_{\mathrm{\delta }}}$ of vertical determination are shown. An anomalous observation signal $Z_V\left(t\right)$ was simulated. In the case of a sequential EKF (20)–(26), estimation errors are shown in Figure 5, and for a robust filter (97)–(98) and (100)–(117), such errors are presented in Figure 6.

Figure 5.

Estimates obtained by the EKF.

Figure 6.

Estimates obtained by the robust filter.

It can be seen that in the absence of protection from anomalous observations, the divergence of the EKF is possible when the actual estimation errors $\mathrm{\Delta \delta }$ considerably differ from their predicted mean square values ${\mathrm{\sigma }}_{\mathrm{\delta }}$.

9. Conclusions

In this chapter, some approaches to monitoring and diagnosis of dynamic systems using the extended Kalman filter were presented. The considered approaches are based on the application of the DS error model, the formation and processing of residuals between the observed parameters and their predicted values. The statistical properties of residuals make it possible to form diagnostic parameters and the corresponding criteria of agreement. For the considered criteria, formalized tolerances can be generated that make it possible to detect failures and failures in the DS.

The conducted studies have confirmed the following possibilities of the proposed approaches.

• Detection and exclusion or processing with certain confidence coefficients of anomalous observations.

• Selection of outliers against the background of gradual biases of observations.

• Improving the reliability of monitoring based on combined goodness-of-fit tests.

• Detection of faults in observable dynamical systems with a depth of the component of the state vector on a basis of the joint procedures of optimal filtering and smoothing of experimental data.