1 Introduction

During the last three decades, the estimation problem of input noise has become an active field in industry and has a wide range of applications in fault detection, petroleum prospecting, image restoration, speech processing, and so forth [16]. The Kalman deconvolution filter approach applied to the reflection coefficient sequence estimate in oil exploration was first proposed in [3]. Later on, the optimal frequency domain deconvolution estimator was derived by employing the polynomial method [7]. Based on the innovation analysis approach, a unified white noise estimator design approach for the autoregressive moving average innovation model was given in [8]. Recently, deconvolution studies have mainly focused on multisensor systems [9, 10] and packet dropout systems [1113]. Because of the limitation of the communication bandwidth, service capacity, and carrying capacity of the network control systems, packet dropouts inevitably exist in the data transmission, which lead to the performance degradation or even instability of the control systems [14]. If it cannot be detected and processed in a timely manner, then it inevitably results in huge losses of personnel and property. Besides, it is worth noting that the input Gaussian noise hypothesis is the basis of the results mentioned, but the non-Gaussian input noise is widespread in numerous important applications [1519]. Inspired by this discussion, we design a novel nonlinear input noise observer for linear stochastic discrete systems with packet dropouts.

The estimation theory of non-Gaussian systems has become more and more influential in many technical fields, and therefore numerous meaningful attempts have been devoted to it [15, 2023]. The non-Gaussian input noise and system state joint estimator was presented for discrete-time nonlinear non-Gaussian systems in [20], where the state posterior distribution was iteratively computed by utilizing the Gaussian sum filtering, and the noise parameter posterior distribution was calculated by applying the variational Bayesian method, respectively. By selecting a suitable system noise a particle estimator with higher efficiency was constructed for nonlinear non-Gaussian systems in [21]. Based on the mixed \(l_{1}\) and \(l_{2}\) norm minimum performance index, an adaptive filter was presented for a class of systems in non-Gaussian and nonlinear manner in [22], where the tuning factor γ was determined by using the projection statistics algorithm. A Tobit Kalman-like estimator was proposed by converting the system with time-correlated and non-Gaussian noises into a Gaussian system with unknown noises covariances in [23]. In [24], by exploiting the Kronecker power of the system states and measurements a nonlinear estimator (called polynomial filtering) was first given for linear non-Gaussian systems, and it was proved that the accuracy of the nonlinear algorithm was higher than that of the classical Kalman filter. Later on, the result in [24] was successfully applied to time-varying systems [25], stochastic bilinear systems [26], general nonlinear stochastic systems [27], and multisensor systems with uncertain observations [28]. In this paper, we use the polynomial filtering theory to investigate the second-order polynomial estimator design problem for a class of packet dropout systems in non-Gaussian manner.

Until now, the research of estimator design of systems in packet dropouts and non-Gaussian manner mainly focuses on state estimate case, but the input noise estimate case is rarely reported. In this study, we propose a design method of input noise nonlinear estimator for this kind of systems. We formulate a new augmented system based on the given system with stochastic packet dropouts by utilizing the original and second-order Kronecker products of the states and observations in the original system. Moreover, we derive the stochastic characteristics of the augmented system by employing some Kronecker algebra rules. Then we adopt the classical Kalman projection theory to produce the recursive second-order polynomial non-Gaussian noise estimator. The main innovations and characteristics of this paper are as follows: (i) to the best of our knowledge, the nonlinear (second-order) estimation of input noise for linear stochastic systems with packet dropouts is investigated for the first time; (ii) a recursive analytical solution of the nonlinear (second-order) estimation of input noise is presented, which can not only realize real-time noise signal estimation, but also has more theoretical significance than the numerical solution.

The rest of our work is arranged as follows. In Sect. 2, we introduce the linear, non-Gaussian, discrete-time, and packet dropout state-space model and define the input noise second-order minimum variance estimation issue. In Sect. 3, we obtain the second-order least-squares filter and fixed-lag smoother of the input noise based on the parameters and stochastic characteristics of a fictitious second-order state-space model. In Sect. 4, we propose a linear input noise recursive estimator by utilizing the classical innovation orthogonal projection theory and compare the performances of linear and second-order noise estimators using an example. Finally, we summarize the research results.

2 Problem statement and preliminaries

Let us introduce the packet dropout state-space model to be discussed:

$$\begin{aligned} \textstyle\begin{cases} x(i+1)=\mathcal{A}(i)x(i)+\mathcal{B}(i)w(i), \\ y(i)=\lambda (i)\mathcal{H}(i)x(i)+\mathcal{D}(i)w(i), \\ z(i)=\mathcal{L}(i)w(i), \\ x(0)=x_{0}, \end{cases}\displaystyle \end{aligned}$$
(1)

where \(i\in \mathbb{N}\) is the discrete-time index, \(x(i)\in \mathbb{R}^{n}\) is the unknown system state, \(y(i)\in \mathbb{R}^{m}\) is the system measurement, \(z(i)\in \mathbb{R}^{q}\) is the noise combination to be used for estimation, and \(\mathcal{A}(i)\in \mathbb{R}^{n\times n}\), \(\mathcal{H}(i)\in \mathbb{R}^{m\times n}\), \(\mathcal{B}(i)\in \mathbb{R}^{n\times r}\), \(\mathcal{D}(i)\in \mathbb{R}^{m\times r}\), and \(\mathcal{L}(i)\in \mathbb{R}^{q\times r}\) are real known time-varying matrices.

Assumption 1

The initial condition \(x_{0}\in \mathbb{R}^{n}\) is a non-Gaussian process independent of \(\{w(i)\}\) that satisfies

$$\begin{aligned} \textstyle\begin{cases} E\{x_{0}\}=0, \\ E\{x_{0}^{[k]}\}=\varPsi _{0,k}, \quad k=2,3,4. \end{cases}\displaystyle \end{aligned}$$
(2)

Assumption 2

The variable \(w(i)\in \mathbb{R}^{r}\) is a zero-mean real vector-valued finite known fourth-order non-Gaussian process satisfying

$$\begin{aligned} \textstyle\begin{cases} E\{w(i)\}=0, \\ E\{w^{[k]}(i)\}=\varPsi _{w,k}, \quad k=2,3,4. \end{cases}\displaystyle \end{aligned}$$
(3)

According to Assumption 2 and the stack operation, we have

$$\begin{aligned} st^{-1}\varPsi _{w,2}=E\bigl\{ w(i)\cdot w^{T}(i)\bigr\} . \end{aligned}$$
(4)

Assumption 3

The binary (0 or 1) stochastic sequence \(\{\lambda (i)\}\in \mathbb{R}\) obeys the following probability distribution:

$$\begin{aligned} \textstyle\begin{cases} \mathrm{Pr}\{\lambda (i)=1\}=p(i), \\ \mathrm{Pr}\{\lambda (i)=0\}=1-p(i), \end{cases}\displaystyle \end{aligned}$$

where \(p(i)\in \mathbb{R}\) is greater than or equal to zero.

The second-order non-Gaussian noise estimation under investigation is as follows: given \(s\in \mathbb{N}\) and based on a measurement sequence \(\{\{y(j)\}_{j=0}^{i+s}\}\), find a second-order polynomial estimator \(\hat{z}(i|i+s)\) of \(z(i)\) that minimizes the mean-squared estimation error.

Remark 1

Similarly to the Kalman filter case, the second-order polynomial estimator \(\hat{z}(i|i+s)\) is a second-order filter when \(s=0\) and a second-order fixed-order smoother when \(s>0\).

3 Main results

In this section, we first construct a second-order extended stochastic system based on Kronecker algebra and non-Gaussian noise hypothesis. Then we design the second-order non-Gaussian noise filter and smoother by utilizing the stochastic properties of the extended system and projection formula.

3.1 The second-order extended system

First, the second-order state vector and measurement vector are defined as follows:

$$\begin{aligned}& x_{e}(i)\triangleq \begin{bmatrix} x(i) \\ x^{[2]}(i)\end{bmatrix}-E \begin{bmatrix} x(i) \\ x^{[2]}(i)\end{bmatrix}, \end{aligned}$$
(5)
$$\begin{aligned}& y_{e}(i)\triangleq \begin{bmatrix} y(i) \\ y^{[2]}(i)\end{bmatrix}-E \begin{bmatrix} y(i) \\ y^{[2]}(i)\end{bmatrix}. \end{aligned}$$
(6)

By (1), (3), and (4) we arrive at

$$\begin{aligned} x^{[2]}(i+1) =&\bigl[\mathcal{A}(i)x(i)+\mathcal{B}(i)w(i)\bigr] \otimes \bigl[ \mathcal{A}(i)x(i)+\mathcal{B}(i)w(i)\bigr] \\ =&\mathcal{A}^{[2]}(i)x^{[2]}(i)+\bigl(\mathcal{A}(i)x(i) \bigr)\otimes \bigl( \mathcal{B}(i)w(i)\bigr) \\ &{}+\bigl(\mathcal{B}(i)w(i)\bigr)\otimes \bigl(\mathcal{A}(i)x(i)\bigr)+ \mathcal{B}^{[2]}(i)w^{[2]}(i) \\ =&\mathcal{A}^{[2]}(i)x^{[2]}(i)+l(i), \end{aligned}$$
(7)

where

$$\begin{aligned} l(i) =&\bigl(\mathcal{A}(i)x(i)\bigr)\otimes \bigl(\mathcal{B}(i)w(i)\bigr)+ \mathcal{B}^{[2]}(i)w^{[2]}(i)+\bigl( \mathcal{B}(i)w(i)\bigr) \otimes \bigl(\mathcal{A}(i)x(i)\bigr) \\ =&\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl( \mathcal{A}(i)x(i)\bigr)\otimes \bigl( \mathcal{B}(i)w(i)\bigr)\bigr]+ \mathcal{B}^{[2]}(i)w^{[2]}(i) \\ =&\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr) \bigl( \bigl(\mathcal{A}(i)x(i)\bigr)\otimes I\bigr) \bigl( \mathcal{B}(i)w(i)\bigr)+ \mathcal{B}^{[2]}(i)w^{[2]}(i). \end{aligned}$$

To construct a fictitious second-order state-space model, we further calculate the following quadratic measurement output:

$$\begin{aligned} y^{[2]}(i) =&\lambda ^{[2]}(i)\mathcal{H}^{[2]}(i)x^{[2]}(i)+ \bigl( \lambda (i)\mathcal{H}(i)x(i)\bigr)\otimes \bigl(\mathcal{D}(i)w(i)\bigr) \\ &{}+\bigl(\mathcal{D}(i)w(i)\bigr)\otimes \bigl(\lambda (i)\mathcal{H}(i)x(i) \bigr)+ \mathcal{D}^{[2]}(i)w^{[2]}(i) \\ =&\lambda ^{[2]}(i)\mathcal{H}^{[2]}(i)x^{[2]}(i)+ \bigl(\mathcal{I}+ \mathcal{O}_{m,m}^{T}\bigr)\bigl[\bigl( \lambda (i)\mathcal{H}(i)x(i)\bigr)\otimes \bigl( \mathcal{D}(i)w(i)\bigr)\bigr] \\ &{}+\mathcal{D}^{[2]}(i)w^{[2]}(i). \end{aligned}$$
(8)

The second-order extended system according to the following equations is yielded by substituting (1) and (7) into (5) and substituting (1) and (8) into (6):

$$\begin{aligned} \textstyle\begin{cases} x_{e}(i+1)=\mathcal{A}_{e}(i)x_{e}(i)+u_{e}(i), \\ y_{e}(i)=\mathcal{D}_{e}(i)\mathcal{H}_{e}(i)x_{e}(i)+w_{e}(i), \end{cases}\displaystyle \end{aligned}$$
(9)

where the system matrix parameters are

$$\begin{aligned} \textstyle\begin{cases} \mathcal{A}_{e}(i)=\operatorname{diag} \{\mathcal{A}(i),\mathcal{A}^{[2]}(i)\}, \\ \mathcal{D}_{e}(i)=\operatorname{diag} \{\lambda (i)\mathcal{I},\lambda ^{[2]}(i) \mathcal{I}\}, \\ \mathcal{H}_{e}(i)=\operatorname{diag} \{\mathcal{H}(i),\mathcal{H}^{[2]}(i)\}, \end{cases}\displaystyle \end{aligned}$$

and

{ u e ( i ) = [ u e , 11 ( i ) u e , 21 ( i ) ] , u e , 11 ( i ) = B ( i ) w ( i ) , u e , 21 ( i ) = ( I + O n , n T ) [ ( A ( i ) x ( i ) ) ( B ( i ) w ( i ) ) ] + B [ 2 ] ( i ) ( w [ 2 ] ( i ) Ψ N , 2 ) ,
(10)
{ w e ( i ) = [ w e , 11 ( i ) w e , 21 ( i ) ] , w e , 11 ( i ) = D ( i ) w ( i ) + [ λ ( i ) p ( i ) ] H ( i ) E [ x ( i ) ] , w e , 21 ( i ) = ( I + O m , m T ) [ ( λ ( i ) H ( i ) x ( i ) ) ( D ( i ) w ( i ) ) ] w e , 21 ( i ) = + [ λ [ 2 ] ( i ) p ( i ) ] H [ 2 ] ( i ) E [ x [ 2 ] ( i ) ] w e , 21 ( i ) = + D [ 2 ] ( i ) ( w [ 2 ] ( i ) Ψ w , 2 ) .
(11)

Then, from the preceding results and Kronecker algebra we come to the following conclusion.

Lemma 1

The driving noises \(\{u_{e}(i)\}\)and measurement noises \(\{w_{e}(i)\}\)in (9), \({i\in \mathbb{N}}\), satisfy the following conditions:

$$\begin{aligned}& E\bigl\{ u_{e}(i)\bigr\} =E\bigl\{ w_{e}(i)\bigr\} =0, \\& E\bigl\{ u_{e}(i)\cdot {u_{e}^{T}(j)}\bigr\} = \mathcal{Q}(i)\delta _{ij}= \begin{bmatrix} \mathcal{Q}_{11}(i)&\mathcal{Q}_{12}(i) \\ \mathcal{Q}_{12}^{T}(i)&\mathcal{Q}_{22}(i)\end{bmatrix}\delta _{ij}, \end{aligned}$$
(12)
$$\begin{aligned}& E\bigl\{ w_{e}(i)\cdot {w_{e}^{T}(j)}\bigr\} = \mathcal{R}(i)\delta _{ij}= \begin{bmatrix} \mathcal{R}_{11}(i)&\mathcal{R}_{12}(i) \\ \mathcal{R}_{12}^{T}(i)&\mathcal{R}_{22}(i)\end{bmatrix}\delta _{ij}, \end{aligned}$$
(13)
$$\begin{aligned}& E\bigl\{ u_{e}(i)\cdot {w_{e}^{T}(j)}\bigr\} = \mathcal{S}(i)\delta _{ij}= \begin{bmatrix} \mathcal{S}_{11}(i)&\mathcal{S}_{12}(i) \\ \mathcal{S}_{21}(i)&\mathcal{S}_{22}(i)\end{bmatrix}\delta _{ij}, \end{aligned}$$
(14)

where the block matrices in (12)(14) are

$$\begin{aligned}& \mathcal{Q}_{11}(i)=\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{B}^{T}(i), \\& \mathcal{Q}_{12}(i)=st^{-1}\bigl[\mathcal{B}^{[3]}(i) \varPsi _{w,3}\bigr], \\& \begin{aligned} \mathcal{Q}_{22}(i)&=\bigl(\mathcal{I}+ \mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl( \mathcal{A}(i)\varPi _{0}(i)\mathcal{A}^{T}(i)\bigr)\otimes \bigl(\mathcal{B}(i) \bigl(st^{-1} \varPsi _{w,2}\bigr)\mathcal{B}^{T}(i) \bigr)\bigr] \\ &\quad {}\times \bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T} \bigr)+\mathcal{B}^{[2]}(i)\bigl\{ st^{-1}\bigl( \varPsi _{w,4}-\varPsi _{w,2}^{[2]}\bigr)\bigr\} { \mathcal{B}^{[2]}}^{T}(i), \end{aligned} \\& \mathcal{R}_{11}(i)=\mathcal{D}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i), \end{aligned}$$
(15)
$$\begin{aligned}& \mathcal{R}_{12}(i)=st^{-1}\bigl[\mathcal{D}^{[3]}(i) \varPsi _{w,3}\bigr], \end{aligned}$$
(16)
$$\begin{aligned}& \begin{aligned}[b] \mathcal{R}_{22}(i)&=\bigl(\mathcal{I}+ \mathcal{O}_{m,m}^{T}\bigr)\bigl[\bigl(p(i) \mathcal{H}(i)\varPi _{0}(i)\mathcal{H}^{T}(i)\bigr) \\ &\quad {}\otimes \bigl(\mathcal{D}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i)\bigr)\bigr]\bigl( \mathcal{I}+ \mathcal{O}_{m,m}^{T}\bigr) \\ &\quad {}+\mathcal{D}^{[2]}(i)\bigl\{ st^{-1}\bigl(\varPsi _{w,4}-\varPsi _{w,2}^{[2]}\bigr)\bigr\} { \mathcal{D}^{[2]}}^{T}(i) \\ &\quad {}+\bigl(p(i)-p^{2}(i)\bigr)\mathcal{H}^{[2]}(i)E \bigl[x^{[2]}(i)\bigr]E^{T}\bigl[x^{[2]}(i)\bigr]{ \mathcal{H}^{[2]}}^{T}(i), \end{aligned} \end{aligned}$$
(17)
$$\begin{aligned}& \mathcal{S}_{11}(i) =\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i), \end{aligned}$$
(18)
$$\begin{aligned}& \mathcal{S}_{12}(i) =st^{-1}\bigl\{ \bigl[ \mathcal{D}^{[2]}(i)\otimes { \mathcal{B}(i)}\bigr]\varPsi _{w,3} \bigr\} , \end{aligned}$$
(19)
$$\begin{aligned}& \mathcal{S}_{21}(i) =st^{-1}\bigl\{ \bigl[\mathcal{D}(i) \otimes {\mathcal{B}^{[2]}(i)}\bigr] \varPsi _{w,3}\bigr\} , \end{aligned}$$
(20)
$$\begin{aligned}& \begin{aligned}[b] \mathcal{S}_{22}(i) &=\bigl(\mathcal{I}+ \mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl(p(i) \mathcal{A}(i)\varPi _{0}(i)\mathcal{H}^{T}(i)\bigr) \\ &\quad {}\otimes \bigl(\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i)\bigr)\bigr]\bigl( \mathcal{I}+ \mathcal{O}_{m,m}^{T}\bigr) \\ &\quad {}+\mathcal{B}^{[2]}(i)\bigl\{ st^{-1}\bigl(\varPsi _{w,4}-\varPsi _{w,2}^{[2]}\bigr)\bigr\} { \mathcal{D}^{[2]}}^{T}(i), \end{aligned} \\& \varPi _{0}(i+1) =\mathcal{A}(i)\varPi _{0}(i) \mathcal{A}^{T}(i)+ \mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{B}^{T}(i), \\& E\bigl[x^{[2]}(i+1)\bigr] =\mathcal{A}^{[2]}(i)E \bigl[x^{[2]}(i)\bigr]+\mathcal{B}^{[2]}(i) \varPsi _{w,2}, \\& E\bigl[x^{[2]}(0)\bigr] =\varPsi _{x,2}. \end{aligned}$$
(21)

Proof

From the hypotheses of state-space model (1) it is known that the additive process \(w(i)\) is independent of \(x_{0}\) and satisfies

$$\begin{aligned} E\bigl\{ w(i)\bigr\} =0, \qquad E\bigl\{ w^{[k]}(i)\bigr\} =\varPsi _{w,k}, \quad k=2,3,4. \end{aligned}$$

From (10) and (11) we can derive the following formulas:

$$\begin{aligned}& E\bigl\{ u_{e}(i)\bigr\} =E\bigl\{ w_{e}(i)\bigr\} =0, \quad {i\in \mathbb{N}}, \\& E\bigl\{ u_{e}(i)\cdot {u_{e}^{T}(j)}\bigr\} =0, \quad {i\neq {j}}, \\& E\bigl\{ w_{e}(i)\cdot {w_{e}^{T}(j)}\bigr\} =0, \quad {i\neq {j}}. \end{aligned}$$

Combining (10) and (12), \(\mathcal{Q}_{11}(i)\) can be given by

$$\begin{aligned} \mathcal{Q}_{11}(i) =&E \bigl\{ \mathcal{B}(i)w(i)w^{T}(i) \mathcal{B}^{T}(i) \bigr\} \\ =&\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr) \mathcal{B}^{T}(i). \end{aligned}$$

Substituting (10) into \(\mathcal{Q}(i)\) and utilizing \(E\{w(i)\cdot x_{0}^{T}\}=0\), we obtain

$$\begin{aligned} \mathcal{Q}_{12}(i) =&E\bigl\{ \mathcal{B}(i)w(i)\bigl[ \mathcal{B}^{[2]}(i) \bigl(w^{[2]}(i)- \varPsi _{w,2} \bigr)\bigr]^{T}\bigr\} \\ =&E \bigl\{ st^{-1}\bigl\{ \bigl[\mathcal{B}^{[2]}(i) \bigl(w^{[2]}(i)-\varPsi _{w,2}\bigr)\bigr] \otimes \bigl( \mathcal{B}(i)w(i)\bigr)\bigr\} \bigr\} \\ =&E \bigl\{ st^{-1}\bigl\{ \bigl[\mathcal{B}^{[2]}(i)\otimes { \mathcal{B}(i)}\bigr] \bigl[\bigl(w^{[2]}(i)- \varPsi _{w,2}\bigr) \otimes {w(i)}\bigr]\bigr\} \bigr\} \\ =&st^{-1}\bigl[\mathcal{B}^{[3]}(i)\varPsi _{w,3} \bigr]. \end{aligned}$$
(22)

Moreover, we have

$$\begin{aligned} \mathcal{Q}_{22}(i) =&E \bigl\{ \bigl[\bigl(\mathcal{I}+ \mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl( \mathcal{A}(i)x(i) \bigr)\otimes \bigl(\mathcal{B}(i)w(i)\bigr)\bigr] \\ &{}+\mathcal{B}^{[2]}(i) \bigl(w^{[2]}(i)-\varPsi _{w,2}\bigr) \bigr] \bigl[\bigl(\mathcal{I}+ \mathcal{O}_{n,n}^{T} \bigr)\bigl[\bigl(\mathcal{A}(i)x(i)\bigr)\otimes \bigl(\mathcal{B}(i)w(i)\bigr) \bigr] \\ &{}+\mathcal{B}^{[2]}(i) \bigl(w^{[2]}(i)-\varPsi _{w,2}\bigr) \bigr]^{T} \bigr\} \\ =&E \bigl\{ \bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr) \bigl[\bigl(\mathcal{A}(i)x(i)\bigr) \otimes \bigl(\mathcal{B}(i)w(i)\bigr)\bigr] \\ &{}\times \bigl[\bigl(\mathcal{A}(i)x(i)\bigr)\otimes \bigl(\mathcal{B}(i)w(i) \bigr)\bigr]^{T}\bigl( \mathcal{I}+\mathcal{O}_{n,n}^{T} \bigr)^{T} \bigr\} \\ &{}+E \bigl\{ \mathcal{B}^{[2]}(i) \bigl(w^{[2]}(i)-\varPsi _{w,2}\bigr) \bigl(w^{[2]}(i)- \varPsi _{w,2} \bigr)^{T}{\mathcal{B}^{[2]}}^{T}(i) \bigr\} \\ =&E \bigl\{ \bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr) \bigl[\bigl(\mathcal{A}(i)x(i)x^{T}(i) \mathcal{A}^{T}(i) \bigr)\otimes \bigl(\mathcal{B}(i)w(i)w^{T}(i)\mathcal{B}^{T}(i) \bigr)\bigr] \bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T} \bigr)^{T} \bigr\} \\ &{}+E \bigl\{ \mathcal{B}^{[2]}(i)st^{-1}\bigl[ \bigl(w^{[2]}(i)-\varPsi _{w,2}\bigr) \otimes \bigl(w^{[2]}(i)- \varPsi _{w,2}\bigr)\bigr]{\mathcal{B}^{[2]}}^{T}(i) \bigr\} \\ =&\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl( \mathcal{A}(i)\varPi _{0}(i) \mathcal{A}^{T}(i)\bigr)\otimes \bigl(\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr) \mathcal{B}^{T}(i)\bigr)\bigr]\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T} \bigr) \\ &{}+E \bigl\{ \mathcal{B}^{[2]}(i)st^{-1} \bigl[w^{[4]}(i)-w^{[2]}(i)\otimes \varPsi _{w,2} -\varPsi _{w,2}\otimes {w^{[2]}(i)}+\varPsi _{w,2}^{[2]} \bigr]{\mathcal{B}^{[2]}}^{T}(i) \bigr\} \\ =&\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T}\bigr)\bigl[\bigl( \mathcal{A}(i)\varPi _{0}(i) \mathcal{A}^{T}(i)\bigr)\otimes \bigl(\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr) \mathcal{B}^{T}(i)\bigr)\bigr]\bigl(\mathcal{I}+\mathcal{O}_{n,n}^{T} \bigr) \\ &{}+\mathcal{B}^{[2]}(i)\bigl\{ st^{-1}\bigl(\varPsi _{w,4}-\varPsi _{w,2}^{[2]}\bigr)\bigr\} { \mathcal{B}^{[2]}}^{T}(i). \end{aligned}$$

Employing the similar lines and properties of non-Gaussian noise \(w(i)\), \(x_{0}\), and \(\lambda (i)\), the covariance matrices \(\mathcal{R}_{11}(i)\), \(\mathcal{R}_{12}(i)\), \(\mathcal{R}_{22}(i)\), \(\mathcal{S}_{11}(i)\), \(\mathcal{S}_{12}(i)\), \(\mathcal{S}_{21}(i)\), and \(\mathcal{S}_{22}(i)\) are directly given by (15)–(21). This completes the proof. □

3.2 The input noise second-order least-squares estimator design

First, we define the innovation \(v_{e}(i)\) and its variance matrix \(R_{v_{e}}(i)\) by

$$ \begin{gathered} v_{e}(i)\triangleq y_{e}(i)-p(i)\mathcal{H}_{e}(i) \hat{x}_{e}(i|i-1), \\ R_{v_{e}}(i)\triangleq E\bigl\{ v_{e}(i)v_{e}^{T}(i) \bigr\} . \end{gathered} $$
(23)

Let \(\hat{x}_{e}(i|i-1)\) denote the least mean-square filter of \(x_{e}(i)\) based on the measurements up to \(i-1\). Then time update is presented as follows:

$$\begin{aligned} \hat{x}_{e}(i+1|i) =&\sum_{j=0}^{i}E \bigl\{ x_{e}(i+1)v_{e}^{T}(j)\bigr\} \mathcal{R}_{v_{e}}^{-1}(j)v_{e}(j) \\ =&\mathcal{A}_{e}(i)\hat{x}_{e}(i|i-1)+ \mathcal{K}(i)v_{e}(i), \end{aligned}$$
(24)

where

$$\begin{aligned}& \mathcal{K}(i) =\bigl[p(i)\mathcal{A}_{e}(i)\mathcal{P}(i) \mathcal{H}_{e}^{T}(i)+ \mathcal{S}(i)\bigr] \mathcal{R}_{v_{e}}^{-1}(i), \end{aligned}$$
(25)
$$\begin{aligned}& \mathcal{R}_{v_{e}}(i) =p(i) \bigl(1-p(i)\bigr)\mathcal{H}_{e}(i) \varPi (i) \mathcal{H}_{e}^{T}(i)+\mathcal{R}(i)+p^{2}(i) \mathcal{H}_{e}(i) \mathcal{P}(i)\mathcal{H}_{e}^{T}(i), \end{aligned}$$
(26)

and

$$\begin{aligned}& \begin{aligned} \varPi (i+1) &=E\bigl\{ x_{e}(i+1)x_{e}^{T}(i+1) \bigr\} \\ &=\mathcal{A}_{e}(i)\varPi (i)\mathcal{A}_{e}^{T}(i)+Q(i), \end{aligned} \\& \begin{aligned} \mathcal{P}(i+1) &=E\bigl\{ \bigl(x_{e}(i+1)-\hat{x}_{e}(i+1|i) \bigr) \bigl(x_{e}(i+1)- \hat{x}_{e}(i+1|i) \bigr)^{T}\bigr\} \\ &=\mathcal{A}_{e}(i)\mathcal{P}(i)\mathcal{A}_{e}^{T}(i)+ \mathcal{Q}(i)- \mathcal{K}(i)\mathcal{R}_{v_{e}}(i)\mathcal{K}^{T}(i). \end{aligned} \end{aligned}$$

In addition, the initial conditions of \(\hat{x}_{e}(i|i-1)\) and \(\mathcal{P}(i)\) are given by

$$\begin{aligned}& \hat{x}_{e}(0|-1) = E\bigl\{ x_{e}(0)\bigr\} , \\& \varPi (0) = E\bigl\{ \bigl(x_{e}(0)-E \bigl(x_{e}(0) \bigr) \bigr) \bigl(x_{e}(0)-E \bigl(x_{e}(0) \bigr) \bigr)^{T}\bigr\} , \end{aligned}$$

which can be easily computed by using (2) in Assumption 1.

Then we obtain a second-order non-Gaussian noise deconvolution filter.

Theorem 1

Under the stochastic packet dropout system (1) with Assumptions 13, we propose a filtering estimation of input noise by the following formula:

$$\begin{aligned} \hat{z}(i|i)=\mathcal{T}(i)v_{e}(i), \end{aligned}$$

where \(v_{e}(i)\)is defined by (23), and

{ T ( i ) = [ T 11 ( i ) T 12 ( i ) ] R v e 1 ( i ) , T 11 ( i ) = L ( i ) ( s t 1 Ψ w , 2 ) G T ( i ) , T 12 ( i ) = s t 1 ( ( G [ 2 ] ( i ) L ( i ) ) Ψ w , 3 ,
(27)

and the innovation variance matrix \(\mathcal{R}_{v_{e}}(i)\)is computed by (26).

Proof

It follows from Kalman projection formula that the deconvolution filtering of \(z(i)\) given measurements from time 0 to time i can be calculated by

$$\begin{aligned} \hat{z}(i|i) =&\sum_{j=0}^{i}E\bigl\{ z(i)v_{e}^{T}(j)\bigr\} \mathcal{R}_{v_{e}}^{-1}(j)v_{e}(j) \\ =&E\bigl\{ z(i)v_{e}^{T}(i)\bigr\} \mathcal{R}_{v_{e}}^{-1}(i)v_{e}(i) \\ =&E \bigl\{ \bigl(\mathcal{L}(i)w(i) \bigr) \bigl(w_{e}(i) \bigr)^{T} \bigr\} \mathcal{R}_{v_{e}}^{-1}(i)v_{e}(i) \\ =& \begin{bmatrix} \mathcal{L}(i)(st^{-1}\varPsi _{w,2})\mathcal{G}^{T}(i)&E ( \mathcal{L}(i)w(i)w_{e,21}^{T}(i) )\end{bmatrix}\mathcal{R}_{v_{e}}^{-1}(i)v_{e}(i), \end{aligned}$$

where the derivation of \(E (\mathcal{L}(i)w(i)w_{e,21}^{T}(i) )\) is similar to the calculation process of (22). The proof is complete. □

Before presenting the input noise smoother, we need to define the covariance

$$\begin{aligned} \mathcal{R}_{i+j}^{i}\triangleq {E\bigl\{ w(i),e_{x}(i+j) \bigr\} }, \quad j=1,2, \ldots ,s, \end{aligned}$$
(28)

where \(e_{x}(i+j)=x_{e}(i+j)-\hat{x}_{e}(i+j|i+j-1)\).

Theorem 2

Under the stochastic packet dropout system (1) with Assumptions 13and given \(s\in \mathbb{N}^{+}\), we propose a smoothing estimation of input noise by the following formula:

$$\begin{aligned} \hat{z}(i|i+s) =\mathcal{T}(i)v_{e}(i)+\mathcal{L}(i)\sum _{j=1}^{s} \bigl(p(i+j)\mathcal{R}_{i+j}^{k} \mathcal{H}_{e}^{T}(i+j)\mathcal{R}_{v_{e}}^{-1}(i+j)v_{e}(i+j) \bigr), \end{aligned}$$
(29)

where \(v_{e}(i+j)\) (\(0\leq {j}\leq {s}\)) is defined by (23), \(\mathcal{T}(i)\)is calculated by (27), \(\mathcal{R}_{v_{e}}(i+j)\) (\(1\leq {j}\leq {s}\)) is given in (26), and \(\mathcal{R}_{i+j}^{i}\) (\(1\leq {j}\leq {s}\)) is recursively calculated by the following matrix formula:

$$\begin{aligned} \mathcal{R}_{i+j+1}^{i} =\mathcal{R}_{i+j}^{i} \bigl[\mathcal{A}_{e}(i+j)-p(i+j) \mathcal{K}(i+j)\mathcal{H}_{e}(i+j) \bigr]^{T}, \quad j=1,2,\ldots ,s-1, \end{aligned}$$
(30)

with

$$\begin{aligned} \mathcal{R}_{i+1}^{i} = & \begin{bmatrix} (st^{-1}\varPsi _{w,2})\mathcal{B}^{T}(i)&st^{-1} ( (\mathcal{B}^{[2]}(i) \otimes \mathcal{I} )\varPsi _{w,3} )\end{bmatrix} \\ &{}- \begin{bmatrix} (st^{-1}\varPsi _{w,2})\mathcal{D}^{T}(i)&st^{-1} ( (\mathcal{D}^{[2]}(i) \otimes \mathcal{I} )\varPsi _{w,3} )\end{bmatrix}\mathcal{K}^{T}(i), \end{aligned}$$

and \(\mathcal{K}(i+j)\) (\(0\leq {j}\leq {s-1}\)) is calculated by (25).

Proof

Since \(z(i)\) is uncorrelated with \(\mathcal{L}\{v_{e}(0),v_{e}(1),\ldots ,v_{e}(i-1)\}\), the calculation of the smoother \(\hat{z}(i|i+s)\) is reduced to

$$\begin{aligned} \hat{z}(i|i+s) =&\sum_{j=0}^{i+s}E\bigl\{ z(i)v_{e}^{T}(j)\bigr\} \mathcal{R}_{v_{e}}^{-1}(j)v_{e}(j) \\ =&\sum_{j=0}^{s}E\bigl\{ z(i)v_{e}^{T}(i+j)\bigr\} \mathcal{R}_{v_{e}}^{-1}(i+j)v_{e}(i+j) \\ =&\mathcal{T}(i)v_{e}(i)+\sum_{j=1}^{s}E \bigl\{ z(i)v_{e}^{T}(i+j)\bigr\} \mathcal{R}_{v_{e}}^{-1}(i+j)v_{e}(i+j). \end{aligned}$$
(31)

Substituting (23) into \(E\{z(i)v_{e}^{T}(i+j)\}\) and using the definition in (28), we get the following expression:

$$\begin{aligned} E\bigl\{ z(i)v_{e}^{T}(i+j)\bigr\} =p(i+j)\mathcal{L}(i) \mathcal{R}_{i+j}^{i} \mathcal{H}_{e}^{T}(i+j). \end{aligned}$$
(32)

Then we search for the expression of \(e_{x}(i+j+1)\). It follows from (9), (23), and (24) that

$$\begin{aligned} e_{x}(i+j+1) =&\bigl[\mathcal{A}_{e}(i+j)-p(i+j) \mathcal{K}(i+j) \mathcal{H}_{e}(i+j)\bigr]e_{x}(i+j) \\ &{}-\mathcal{K}(i+j) \bigl(\mathcal{D}_{e}(i+j)-p(i+j)I\bigr) \mathcal{H}_{e}(i+j) \\ &{}\times x_{e}(i+j)+u_{e}(i+j)-\mathcal{K}(i+j)w_{e}(i+j). \end{aligned}$$
(33)

Taking into account that \(w(i)\), \(u_{e}(i+j)\), and \(w_{e}(i+j)\) are uncorrelated to each other, we can obtain (30) by substituting (33) into (28).

Furthermore, fixed-lag smoother expression (29) follows directly by substituting (32) into (31). This completes the proof. □

Remark 2

With the development of information technology and internet of things, many modern industrial systems are often based on network control, and random packet dropouts are inevitable in network control. The existing estimator design theory of conventional linear systems is not suitable for solving such complex problems. Although a design method of second-order estimator for non-Gaussian noise is proposed in [25], it is only suitable for conventional linear stochastic non-Gaussian systems, and it is only a particular case of this paper (that is, when the probability of the random variable \(\lambda (i)\) in system (1) is \(p(i)=1\), the system in this paper is equivalent to the system in [25]). Therefore the conclusion of this paper is more general than that in [25].

4 Linear input noise estimator and numerical example

In this section, by a numerical simulation we show the validity of the presented algorithm. Note that, in the mean-squared sense, the Kalman-based estimator is the best linear estimator when the input noise is assumed to be non-Gaussian [24]. For comparing the linear input noise estimator with the second-order input noise estimator, we first propose the linear deconvolution filtering of \(z(i)\) by applying the Kalman filtering theory in the following equations.

Define the innovation process

$$\begin{aligned} m(i)\triangleq y(i)-\hat{y}(i|i-1) , \end{aligned}$$

and \(\mathcal{R}_{m}(i)\triangleq E\{m(i)\cdot m^{T}(i)\}\), which can be computed by

$$\begin{aligned} \mathcal{R}_{m}(i)=p(i) \bigl(1-p(i)\bigr)\mathcal{H}(i)\varPi _{0}(i)\mathcal{H}^{T}(i)+ \mathcal{D}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i) +p^{2}(i) \mathcal{H}(i)\mathcal{P}_{0}(i) \mathcal{H}^{T}(i), \end{aligned}$$

where

$$\begin{aligned} \varPi _{0}(i+1) \triangleq &E\bigl\{ x(i+1)x^{T}(i+1)\bigr\} \\ =&\mathcal{A}(i)\varPi _{0}(i)\mathcal{A}^{T}(i)+ \mathcal{B}(i) \bigl(st^{-1} \varPsi _{w,2}\bigr) \mathcal{B}^{T}(i) \end{aligned}$$

and

$$\begin{aligned} \mathcal{P}_{0}(i)\triangleq E\bigl\{ \bigl(x(i)-\hat{x}(i|i-1)\bigr) \bigl(x(i)-\hat{x}(i|i-1)\bigr)^{T} \bigr\} . \end{aligned}$$
(34)

Moreover, by applying the projection formula and innovation process we have

$$\begin{aligned} \hat{x}(i+1|i)=\mathcal{A}(i)\hat{x}(i|i-1)+\mathcal{K}_{0}(i)m(i), \end{aligned}$$
(35)

where

$$\begin{aligned} \mathcal{K}_{0}(i)= \bigl(p(i)\mathcal{A}(i)\mathcal{P}_{0}(i) \mathcal{H}^{T}(i)+\mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{D}^{T}(i) \bigr)\mathcal{R}_{m}^{-1}(i). \end{aligned}$$

By (1), (34), and (35) we obtain

$$\begin{aligned} \mathcal{P}_{0}(i+1)=\mathcal{A}(i)\mathcal{P}_{0}(i) \mathcal{A}^{T}(i)+ \mathcal{B}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{B}^{T}(i)-\mathcal{K}_{0}(i) \mathcal{R}_{m}(i)\mathcal{K}^{T}_{0}(i). \end{aligned}$$

Thus the linear input noise filter \(\hat{z}(i|i)\) of \(z(i)\) is computed as follows:

$$\begin{aligned} \hat{z}(i|i)=\mathcal{L}(i) \bigl(st^{-1}\varPsi _{w,2}\bigr) \mathcal{D}^{T}(i) \mathcal{R}_{m}^{-1}(t)m(i), \end{aligned}$$

and the linear smoother \(\hat{z}(i|i+s)\) (\(s\in \mathbb{N}^{+}\)) of \(z(i)\) is calculated by

$$\begin{aligned} \hat{z}(i|i+s)=\hat{z}(i|i)+\mathcal{L}(i)\sum_{j=1}^{s}p(i+j) \varUpsilon _{i+j}^{i}\mathcal{H}^{T}(i+j) \mathcal{R}_{m}^{-1}(i+j)m(i+j), \end{aligned}$$

where

$$\begin{aligned} \varUpsilon _{i+j+1}^{i}=\varUpsilon _{i+j}^{i} \bigl[\mathcal{A}(i+j)-p(i+j) \mathcal{K}_{0}(i+j)\mathcal{H}(i+j) \bigr]^{T}, \quad j=1,2,\ldots ,s-1, \end{aligned}$$

with

$$\begin{aligned} \varUpsilon _{i+1}^{i}=\bigl(st^{-1}\varPsi _{w,2}\bigr)\mathcal{B}^{T}(i)-\bigl(st^{-1} \varPsi _{w,2}\bigr)\mathcal{D}^{T}(i)\mathcal{K}_{0}^{T}(i). \end{aligned}$$

Furthermore, we assume that the system matrix parameters in (1) are as follows:

$$\begin{aligned}& \mathcal{A}(i)=0.7+0.2e^{-0.3i}, \qquad \mathcal{B}(i)=0.2\cos (0.1i), \qquad \mathcal{H}(i)=0.5+0.2\cos (5i), \\& \mathcal{D}(i)=0.4+0.1e^{-7i}, \qquad \mathcal{L}(i)=0.5, \end{aligned}$$

where the system noise vector \(w(i)\) and the initial condition vector \(x_{0}\) are uncorrelated zero-mean non-Gaussian processes with distribution laws given in Table 1.

Table 1 Distribution laws of \(w(i)\) and \(x_{0}\)

The Bernoulli stochastic variable \(\lambda (i)\) is assumed to have a nonzero probability \(p(i)=0.9\), and its variation law is shown in Fig. 1. With the same assumptions, the linear and second-order non-Gaussian noise estimates have been compared in this numerical simulation. The performances of linear estimators and second-order estimators have been plotted in Figs. 25. Obviously, the simulation of this example shows that the second-order estimators are better than the linear estimators.

Figure 1
figure 1

Bernoulli stochastic variable \(\lambda (i)\) when \({p(i)=0.9}\)

Figure 2
figure 2

Signal, linear filter, and second-order filter

Figure 3
figure 3

Linear filter error and second-order filter error

Figure 4
figure 4

Signal, linear four-step lag smoother, and second-order four-step lag smoother

Figure 5
figure 5

Linear four-step lag smoother error and second-order four-step lag smoother error

5 Conclusions

A novel input noise nonlinear (second-order) estimate algorithm is put forward for non-Gaussian stochastic systems with packet dropouts, where the packet dropout characteristic is modeled as a multiplicative binary (0 or 1) distributed stochastic white sequence. By defining the original and second-order Kronecker products of the states and measurements in the original system we introduce a new augmented system. On this basis, we adopt the Kronecker algebra rules to analyze the stochastic characteristics of the augmented system. Then the input noise nonlinear (second-order) estimators are presented in the form of filtering and smoothing. In addition, we emphasize the effectiveness of the proposed algorithm by a numerical simulation.

As a matter of fact, only some basic results have been achieved in the study of complex systems, such as the stability analysis of switched impulsive systems [29] and stochastic stability analysis of nonlinear second-order stochastic differential delay systems [30]. Therefore, it would be of interest to extend the proposed method to investigate more complex switched impulsive systems and nonlinear stochastic delay systems. Another interesting open topic is the \(H_{\infty }\) nonlinear (second-order or higher-order) estimator design problem for systems with packet dropouts.