Recursive Input Reconstruction via Two Kalman Filters and Its Applications

1. Introduction

Inverse problems are important problems in both science and engineering [1,2], such as virtual sensing [3], image processing (e.g., image restoration [4]), sensor linearization [5], and digital predistortion for radio frequency communications [5]. In the point of view of system identification, there are two kinds of inverse problems [6]:

(i)

Reconstruct the system inputs based on the system outputs and the inverse system model, which is also called the input reconstruction problem or the inverse system identification problem, see Figure 1(a).

(ii)

Identify the forward system model based on the input-output data, which is a normal system identification problem, see Figure 1(b).

In this paper, the first type of inverse problems (i.e., the input reconstruction problem) is investigated. Generally, there are two types of approaches for solving the first kind of inverse problems:

(i)

The first approach is to make a direct inversion of the nominal system firstly, and then input reconstruction can be conducted. Denote the transfer function of a discrete-time model as $\left(z\right)$ of which a state-space realization is $(,,,)$, when the inverse of the feedthrough term does not exist, the direct inversion of the model $\left(z\right)$ cannot be conducted [7,8]. In addition, if there exist nonminimum-phase zeros in $\left(z\right)$, an unstable inversion solution will be obtained [9]. So in practical applications of direct inversion approaches,

(ii)

The second approach is to obtain an inverse system model of the nominal system model indirectly, and input reconstruction can then be realized. However, in order to obtain a stable inversion, there are a number of drawbacks in existing approaches:

(a)

For exact stable inversion approaches, an infinite preview of the desired output is needed [10,11,12,13], which is not applicable in practice.

(b)

For the approaches called pseudo-inversion, even though they do not need preview, however, they will encounter other problems such as the difficulty of choosing a suitable basis function [14,15].

(c)

For stable and unstable poles separation-based inversion methods, a finite preview is necessary [16,17].

(d)

For norm-based inversion methods, some of them suffer the problems such as input noise [5] and non-convex optimization [18,19] which often occur in system identification, while for $H_{\infty }$ norm-based methods [9,20,21,22,23], a finite preview is needed.

(e)

For signal modelling-based inversion methods, the input signal, which is to be reconstructed, must be a periodic signal under stationary operating conditions [24,25].

As seen, for some indirect system inversion approaches, even though a stale inversion can be obtained, an infinite or a finite pre-actuation is still needed, which cannot be applied well in practice because sometimes the desired output in unknown. In order to solve the input reconstruction problem in a better way, in this paper an alternative approach is proposed. The presented approach can guarantee the stability of the input reconstructor, and simultaneously the proposed approach does not need any pre-actuation; Moreover, the approach can be applied to stable or unstable, proper or improper systems1 with input to be reconstructed, and there is also no requirement for the type of input and output signals; Furthermore, it does not suffer non-convex or input noise problems.

The remainder of the paper is organized as follows. In Section 2, the modeling of signals with finite-length is introduced, based on which in Section 3 an alternative recursive Kalman filter-based input reconstruction approach is proposed. In Section 4, the performance of the proposed approach is verified and analyzed, and finally conclusions and future perspectives are given in Section 5.

2. Limited-Length Signal Modeling

In this section, the modeling of limited-length signals is illustrated. With the idea of Limited-length signal modeling, an alternative input reconstruction approach is proposed in Section 3.

Given a discrete-time signal $u\left(k\right)\in \mathbb{R}$ with length N and the sampling period $T_{\mathrm{s}}$ in seconds, now assume that the limited-length signal $u\left(k\right)$ is a whole period of a periodic signal $u_{\mathrm{p}}\left(k\right)\in \mathbb{R}$, so within the length N, the signal $u\left(k\right)$ can be represented as the output of the following state-space model [24], i.e.,

$$\{\begin{array}{c}\hfill x_{\mathrm{u}}(k+1)=A_{\mathrm{u}}{x}_{\mathrm{u}}\left(k\right),\\ \hfill u\left(k\right)=C_{\mathrm{u}}{x}_{\mathrm{u}}\left(k\right)+v_{\mathrm{u}}\left(k\right),\end{array}$$

(1)

where $x_{\mathrm{u}}\left(k\right)\in {\mathbb{R}}^{n_{\mathrm{u}}}$ denotes the state vector, the term $v_{\mathrm{u}}\left(k\right)$ denotes the modeling error which is induced by the limited dimension of the matrix $A_{\mathrm{u}}$, and the matrices $A_{\mathrm{u}}$ and $C_{\mathrm{u}}$ can be represented as

$$A_{\mathrm{u}}=\left(\begin{array}{cccc}1& 0& \cdots & 0\\ 0& A_1& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & 0& A_{n_{\mathrm{u}}}\end{array}\right)$$

(2)

and

$$C_{\mathrm{u}}=\left(\begin{array}{cccc}1& C_1& \cdots & C_{n_{\mathrm{u}}}\end{array}\right),$$

(3)

respectively, where denotes a zero vector or a zero matrix, and the individual block entries in these block matrices are

$$A_i=\left(\begin{array}{cc}cos\left(2\pi i{f}_{\mathrm{b}}{T}_{\mathrm{s}}\right)& sin\left(2\pi i{f}_{\mathrm{b}}{T}_{\mathrm{s}}\right)\\ -sin\left(2\pi i{f}_{\mathrm{b}}{T}_{\mathrm{s}}\right)& cos\left(2\pi i{f}_{\mathrm{b}}{T}_{\mathrm{s}}\right)\end{array}\right)$$

(4)

and

$$C_i=\left(\begin{array}{cc}1& 0\end{array}\right),$$

(5)

where $f_{\mathrm{b}}=\frac{1}{N{T}_{\mathrm{s}}}$, for $i=1,2,\dots ,n_{\mathrm{u}}$.

In the end of this section, the following assumption is made.

Assumption A1.

The sequence $\left\{v_{\mathrm{u}}\left(k\right)\right\}$ is assumed to be a white-noise sequence. The covariance function of the white-noise sequence $\left\{v_{\mathrm{u}}\left(k\right)\right\}$ is ${\sigma }^2{\delta }_{kj}$ where ${\sigma }^2$ is a positive and constant scalar value, and ${\delta }_{kj}$ is the Kronecker Delta function depending on two integral numbers k and j [26]:

$${\delta }_{kj}=\{\begin{array}{c}\hfill 0,\mathrm{if}k\ne j,\\ \hfill 1,\mathrm{if}k=j.\end{array}$$

(6)

3. Input Reconstruction Approach

In this section, an approach solving the reconstruction problem of an infinite-length input is derived.

Consider the following discrete-time, linear, time-invariant model which is minimal-realized2 and proper3:

$$\{\begin{array}{c}\hfill A(k+1)=Ax\left(k\right)+Bu\left(k\right),\\ \hfill y\left(k\right)=Cx\left(k\right)+Du\left(k\right)+v\left(k\right),\end{array}$$

(7)

where $u\left(k\right)\in \mathbb{R}$, $y\left(k\right)\in {\mathbb{R}}^p$, and $x\left(k\right)\in {\mathbb{R}}^n$ denote the input, the output, and the state vector, respectively. The matrices A, B, C, D are the state matrix, the input matrix, the output matrix, and the feedthrough matrix, respectively. The term $v\left(k\right)\in {\mathbb{R}}^q$ represents a noise term.

The sampling period of the discrete-time model (2) is the same as the sampling period of the model (1), i.e., $T_{\mathrm{s}}$ in seconds. The model (7) can be stable or unstable.

In the model (7), the input signal $u\left(k\right)$ is with length N, then $u\left(k\right)$ can be modeled, see the modeling process in Section 2. By augmenting the state vector $x_{\mathrm{u}}\left(k\right)$ of the model (1) with the state vector $\left(k\right)$ of the model (7), we can obtain the following augmented model:

$$\{\begin{array}{c}\hfill x_{\mathrm{a}}(k+1)=A_{\mathrm{a}}{x}_{\mathrm{a}}\left(k\right)+v_1\left(k\right),\\ \hfill y\left(k\right)=C_{\mathrm{a}}{x}_{\mathrm{a}}\left(k\right)+v_2\left(k\right),\end{array}$$

(8)

where the state vector

$$x_{\mathrm{a}}\left(k\right)=\left(\begin{array}{c}{x}_{\mathrm{u}}\left(k\right)\\ x\left(k\right)\end{array}\right)\in {\mathbb{R}}^{n_{\mathrm{a}}},$$

(9)

and the matrices

$$A_{\mathrm{a}}=\left(\begin{array}{cc}{A}_{\mathrm{u}}& 0\\ B{C}_{\mathrm{u}}& A\end{array}\right)$$

(10)

and

$$C_{\mathrm{a}}=\left(\begin{array}{cc}D{C}_{\mathrm{u}}& C\end{array}\right),$$

(11)

and the noise terms

$$v_1\left(k\right)=\left(\begin{array}{c}0\\ B{v}_{\mathrm{u}}\left(k\right)\end{array}\right)$$

(12)

and

$$v_2\left(k\right)=D{v}_{\mathrm{u}}\left(k\right)+v\left(k\right).$$

(13)

The following assumptions are made for the above noise terms.

Assumption 2.

The distribution of the initial state vector $x_{\mathrm{a}}\left(0\right)$ is Gaussian.

Assumption A3.

The sequence $\left\{\right(k\left)\right\}$ is assumed to be a white-noise sequence, and the covariance matrix of the white noise sequence $\left\{\right(k\left)\right\}$ is ${\delta }_{kj}$ with a positive and constant matrix.

Assumption A4.

The sequences ${\{}_1\left(k\right)\}$ and ${\{}_2\left(k\right)\}$ are uncorrelated with the initial state vector $x_{\mathrm{a}}\left(0\right)$. The covariance matrices of ${\{}_1\left(k\right)\}$ and ${\{}_2\left(k\right)\}$ are denoted as $Q_{\mathrm{a}}{\delta }_{kj}$ and $R_{\mathrm{a}}{\delta }_{kj}$, respectively. With the values of ${\sigma }^2$ and , the values of $Q_{\mathrm{a}}$ and $R_{\mathrm{a}}$ can be calculated according to (12) and (13).

Below Section 3.1 and Section 3.2 are used to demonstrate the idea of the proposed input reconstruction approach. Section 3.1 illustrates how a limited-length input can be reconstructed, based on which a recursive reconstruction approach of an finite-length input is derived.

3.1. Limited-Length Input Reconstruction

Based on the above analysis and assumptions, the Kalman filter for the model (8) can be implemented [26]. Denote the conceptual time-varying transfer operator of the Kalman filter for the model (8) as $K_{\mathrm{f}}(q^{-1},k)$, and then

$${\widehat{x}}_{\mathrm{a}}\left(k\right)=K_{\mathrm{f}}(q^{-1},k)\left(k\right),$$

(14)

where the notation “$\widehat{ }$” denotes the estimate, and then

$$\widehat{u}\left(k\right)=C_{\mathrm{r}}{\widehat{x}}_{\mathrm{a}}\left(k\right),$$

(15)

where $C_{\mathrm{r}}=\left(\begin{array}{cc}{C}_{\mathrm{u}}& 0\end{array}\right)$, for $k=0,1,\dots ,N$.

According to (9), the dimension of the state vector $x_{\mathrm{a}}\left(k\right)$ is

$$n_{\mathrm{a}}=n+2f_{\mathrm{m}}N{T}_{\mathrm{s}}+1,$$

(16)

where $f_{\mathrm{m}}$ in Hz denotes the frequency corresponding to the largest frequency component number (i.e., $n_{\mathrm{u}}$) of the signal which is denoted as

$$u_{\mathrm{m}}\left(k\right)=u\left(k\right)-v_{\mathrm{u}}\left(k\right).$$

(17)

3.2. Recursive Input Reconstruction Algorithm

In Section 3.1, the Kaman filter for the model (8) can be successfully used for the reconstruction of the input signal $u\left(k\right)$ with length N, and it is clear that the drawbacks in existing indirect system inversion approaches mentioned in Section 1 can be avoided, however, the Kalman filter is merely effective in a finite-time horizon. Below a reconstruction approach of the input signal $u\left(k\right)$ with an infinite length is derived by involving another Kalman filter for the model (8), i.e., based on a collaboration of two Kalman filters, the reconstruction problem of an infinite-length signal can be solved.

The specific mechanism of how to use the two Kalman filters alternatively to solve the reconstruction problem of an infinite-length signal can be illustrated by one timing diagram in Figure 2 and three logic blocks in Figure 3:

(i)

Timing Diagram

In Figure 2, the start and end time points of the two Kalman filters ${\mathit{K}}_1$ and ${\mathit{K}}_2$ are displayed. In more detail, for ${\mathit{K}}_1$, it starts from the step $k=0$ and stops at the step $k=N-1$, and then starts again from $k=N$ and stops at $k=2N-1$, and so forth. While ${\mathit{K}}_2$ starts from $k=c$ and stops at $k=N-1+c$, and then starts again from $k=N+c$ and stops at $k=2N-1+c$, and so on. The two Kalman filters have overlapped working periods (e.g., the shading part in Figure 2), with the value c, the two Kalman filters ${\mathit{K}}_1$ and ${\mathit{K}}_2$ can be implemented alternatively by using the logic blocks shown in Figure 3 such that both transient and finite-length problem can be solved.

(ii)

Logic Block 1

The logic block 1 is used for initializing the prediction process in ${\mathit{K}}_1$, i.e., at steps $k_{1}N$, for $k_1=0,1,2,\dots$, the vector $x_{\mathrm{a}}\left(k\right)$ and the matrix $\left(k\right)$ are forced to be ${\mathbf{0}}_{n_{\mathrm{a}}}$ (i.e., an $n_{\mathrm{a}}$-by-1 zero vector) and ${\mathit{I}}_{n_{\mathrm{a}}}$ (an $n_{\mathrm{a}}$-by-$n_{\mathrm{a}}$ identity matrix) selected from the initial value bank, respectively.

(iii)

Logic Block 2

The logic block 2 is used for the initialization of ${\mathit{I}}_2$, i.e., at steps $k_{2}N+c$, for $k_2=0,1,2,\dots$, the vector $x_{\mathrm{a}}\left(k\right)$ and the matrix $\left(k\right)$ are forced to be ${\mathbf{0}}_{n_{\mathrm{a}}}$ and ${\mathit{I}}_{n_{\mathrm{a}}}$ selected from the initial value bank, respectively.

(iv)

Logic Block 3

The logic block 3 is used for reconstructing the input signal $\widehat{u}\left(k\right)$. As seen in (15), the reconstructed input signal $\widehat{u}\left(k\right)$ can be calculated by using the estimate ${\widehat{u}}_{\mathrm{a}}\left(k\right)$, based on which the specific idea behind the logic block 3 is illustrated in (18):

$$\widehat{u}\left(k\right)=\{\begin{array}{c}\hfill C_{\mathrm{r}}{\widehat{x}}_{\mathrm{a}}^{\left(1\right)}\left(k\right),\mathrm{if}k\in A_1,\\ \hfill C_{\mathrm{r}}{\widehat{x}}_{\mathrm{a}}^{\left(2\right)}\left(k\right),\mathrm{if}k\in A_2,\end{array}$$

(18)

where ${\widehat{x}}_{\mathrm{a}}^{\left(1\right)}\left(k\right)$ denotes the state vector $x_{\mathrm{a}}\left(k\right)$ estimated by ${\mathit{K}}_1$ while ${\widehat{x}}_{\mathrm{a}}^{\left(2\right)}\left(k\right)$ represents the state vector $x_{\mathrm{a}}\left(k\right)$ estimated by ${\mathit{K}}_2$, the sets $A_1$ and $A_2$ are respectively represented as

$$A_1=[0,N-1-d]\cup [k_{3}N+c-d,(k_3+1)N-1-d],$$

(19)

for $k_3=1,2,3,\dots$, and

$$A_2=[k_{4}N-d,k_{4}N+c-1-d],$$

(20)

for $k_4=1,2,3,\dots$, where d is a positive integer, the reason why d is involved is that in practice the part of not interest in $u\left(k\right)$ is usually unknown such that the signal model (1) is not accurately enough to represent the signal $u\left(k\right)$ with length N.

Based on the above analysis, the algorithm of reconstructing the input signal $u\left(k\right)$ with infinite length can be demonstrated in Algorithm Section 3.2, in which ${\mathit{K}}_1(q^{-1},k)$ and ${\mathit{K}}_2(q^{-1},k)$ denote the conceptual time-varying transfer operators of the Kalman filters ${\mathit{K}}_1$ and ${\mathit{K}}_2$, respectively, and the superscripts $\left(1\right)$ and $\left(2\right)$ are used for distinguishing between ${\mathit{K}}_1$ and ${\mathit{K}}_2$.

4. Numerical Simulation and Analysis

In this section, the proposed input reconstruction approach is validated using two examples. In the first example, a randomly generated system is used for validating the effectiveness of the proposed approach. In the second example, a mass-spring-damper system is used to validate the effectiveness

4.1. Input Reconstruction of a Randomly Generated System

Firstly, the system matrices for the model (7) are randomly generated:

$$\mathit{A}=\left(\begin{array}{cccc}-0.4055& -0.5366& 0.1323& 0.2417\\ 0.0904& -0.3823& 0.5593& 0.2120\\ -0.5795& -0.2297& -0.5002& 0.1414\\ -0.1410& 0.1588& -0.2798& -0.2501\end{array}\right),$$

(21)

$$\mathit{B}={\left(\begin{array}{cccc}-1.0164& 2.9668& 0.0826& 0.0870\end{array}\right)}^{\mathrm{T}},$$

(22)

$$\mathit{C}=\left(\begin{array}{cccc}-0.3544& 0.7639& 0.9837& 1.4688\\ -1.3393& 0.4579& 0& 0.1052\end{array}\right),$$

(23)

and

$$\mathit{D}={\left(\begin{array}{cc}0& 0\end{array}\right)}^{\mathrm{T}}.$$

(24)

Secondly, choose a white noise sequence with the covariance function ${\sigma }_{\mathrm{w}}^2{\delta }_{kj}$ with

$${\sigma }_{\mathrm{w}}^2=1\times {10}^{-2},$$

(25)

and then pass it through a 7th-order Butterworth low-pass filter with a cutoff frequency of 300 Hz, afterwards choose the filtered white noise as the input $u\left(k\right)$. Additionally,

Thirdly, set ${\sigma }^2=1\times {10}^{-2}$, and set the value of to be

$$\left(\begin{array}{cc}1\times {10}^{-6}& 0\\ 0& 1\times {10}^{-6}\end{array}\right).$$

Fourthly, choose the sampling period $T_{\mathrm{s}}$ as $1\times {10}^{-4}$ seconds, and then based on the model (7), the simulated output $\left(k\right)$ can be obtained.

Finally, with the simulated output $\left(k\right)$, based on using Algorithm Section 3.2, in which $f_{\mathrm{m}}=300$ Hz, $N=1000$, $c=200$, $d=50$, and $T_{\mathrm{s}}=1\times {10}^{-4}$ seconds, the reconstructed input signal $\widehat{u}\left(k\right)$ can be obtained.

The input reconstruction results from the step 800 to the step 1400 based on using Algorithm Section 3.2 is illustrated in Figure 4. As seen in Figure 4, by comparing the difference between $u\left(k\right)$ and $\widehat{u}\left(k\right)$, it can be concluded that the proposed input reconstruction approach is effective, and it can be used for a signal with infinite length.

In Figure 5, it can be seen that without the logic block 3 in the proposed approach, merely using ${\mathit{K}}_1$ or ${\mathit{K}}_2$ cannot reconstruct the input signal well, especially near the step 1000 and the step 1200.

4.2. Virtual Input Force Sensing

Figure 6 displays a mass-spring-damper system which can be found in many mechanical systems. One common example is the suspension of a car. In this example, the input force u is reconstructed by using the displacement $d_{\mathrm{s}}$ of the mass m, i.e., a virtual input force sensor is formulated using the proposed input reconstruction approach and the data from the displacement sensor measuring the value of $d_{\mathrm{s}}$.

The following numerical values to the variables m, $k_{\mathrm{s}}$, and b in the mass-spring-damper system is shown in Table 1.

According to the above values, the continuous-time transfer function $G\left(s\right)$ between the input u and the output $d_{\mathrm{s}}$ can be obtained:

$$G\left(s\right)=\frac{1}{m{s}^2+bs+k_{\mathrm{s}}}=\frac{1}{s^2+0.2s+1}.$$

(26)

Transform the continuous-time transfer function $G\left(s\right)$ into a discrete-time state-space model with the sampling period $T_s=1\times {10}^{-4}$ seconds, then the system matrices of the state-space model can be obtained:

$$\mathit{A}=\left(\begin{array}{cc}-0.2& -1.0\\ 1& 0\end{array}\right),$$

(27)

$$\mathit{B}={\left(\begin{array}{cc}1& 0\end{array}\right)}^{\mathrm{T}},$$

(28)

$$\mathit{C}=\left(\begin{array}{cc}0& 1\end{array}\right),$$

(29)

and

$$\mathit{D}=0.$$

(30)

Given the input force

$$u\left(k\right)=10sin\left(2\pi f_{1}k{T}_{\mathrm{s}}\right)+5sin\left(2\pi f_{2}k{T}_{\mathrm{s}}\right),$$

(31)

where the frequency $f_1=40$ Hz and the frequency $f_2=60$ Hz.

Set ${\sigma }^2=1\times {10}^{-2}$, and set $=1\times {10}^{-6}$, and choose the sampling period $T_{\mathrm{s}}$ as $1\times {10}^{-4}$ seconds, and then based on the model (7), the simulated output $d_{\mathrm{s}}\left(k\right)$ can be obtained. Finally, with the simulated output $d_{\mathrm{s}}\left(k\right)$, based on using Algorithm Section 3.2, in which $f_{\mathrm{m}}=300$ Hz, $N=1000$, $c=200$, $d=50$, and $T_{\mathrm{s}}=1\times {10}^{-4}$ seconds, the virtual input force sensor output $\widehat{u}\left(k\right)$ can be obtained.

The virtual input force sensing results from the step 800 to the step 1400 based on using Algorithm Section 3.2 is illustrated in Figure 7. As seen in Figure 7, by comparing the difference between $u\left(k\right)$ and $\widehat{u}\left(k\right)$, it can be concluded that the proposed input reconstruction approach is effective in the application of virtual sensing.

In Figure 8, it can be seen that without the logic block 3 in the proposed approach, merely using ${\mathit{K}}_1$ or ${\mathit{K}}_2$ cannot reconstruct the input signal well, especially near the step 1000 (not obvious) and the step 1200.

According to the applications of Algorithm Section 3.2 in the above examples, direct inversion of the model (7) can be avoided, and only the detectability4 of the model (8) should be guaranteed, and in addition, the proposed approach can avoid the drawbacks in existing indirect system inversion approaches mentioned in Section 1.

5. Conclusions and Perspectives

In this paper, an alternative approach for reconstructing infinite-length signal is proposed, and the proposed approach is prospective method for virtual sensing, sensor linearization and so on. The validation results from two examples can illustrate the effectiveness of the proposed approach. According to both theoretical derivation and numerical studies, it can be known that the proposed approach can avoid all the drawbacks existing in current input reconstruction approaches. However, because of the observability problem of the augmented model (8), the proposed approach can be only used for reconstructing a single input or multiple inputs without spectrum overlapping. In the future, the method will be extended to solve the multiple-input reconstruction problem and the problem of input reconstruction of nonlinear systems. Moreover, the effects from the parameters $f_{\mathrm{m}}$, N, c, d, $Q_{\mathrm{a}}$, and $R_{\mathrm{a}}$, and the Kalman filter dimension problem will be explored.