Open AccessArticle

A Two-Filter Approach for State Estimation Utilizing Quantized Output Data

Angel L. Cedeño

^1,2,*

Ricardo Albornoz

^1,2

Rodrigo Carvajal

Boris I. Godoy

⁴

and

Juan C. Agüero

^1,2

Departamento Electrónica, Universidad Técnica Federico Santa María (UTFSM), Av. España 1680, Valparaíso 2390123, Chile

Advanced Center for Electrical and Electronic Engineering, AC3E, Gral. Bari 699, Valparaíso 2390136, Chile

Escuela de Ingeniería Eléctrica, Pontificia Universidad Católica de Valparaíso, Av. Brasil 2147, Valparaíso 2374631, Chile

⁴

Department of Mechanical Engineering, Boston University, Boston, MA 02215, USA

Author to whom correspondence should be addressed.

Sensors 2021, 21(22), 7675; https://doi.org/10.3390/s21227675

Submission received: 27 September 2021 / Revised: 16 November 2021 / Accepted: 16 November 2021 / Published: 18 November 2021

(This article belongs to the Section Sensor Networks)

Download

Browse Figures

Versions Notes

Abstract

Filtering and smoothing algorithms are key tools to develop decision-making strategies and parameter identification techniques in different areas of research, such as economics, financial data analysis, communications, and control systems. These algorithms are used to obtain an estimation of the system state based on the sequentially available noisy measurements of the system output. In a real-world system, the noisy measurements can suffer a significant loss of information due to (among others): (i) a reduced resolution of cost-effective sensors typically used in practice or (ii) a digitalization process for storing or transmitting the measurements through a communication channel using a minimum amount of resources. Thus, obtaining suitable state estimates in this context is essential. In this paper, Gaussian sum filtering and smoothing algorithms are developed in order to deal with noisy measurements that are also subject to quantization. In this approach, the probability mass function of the quantized output given the state is characterized by an integral equation. This integral was approximated by using a Gauss–Legendre quadrature; hence, a model with a Gaussian mixture structure was obtained. This model was used to develop filtering and smoothing algorithms. The benefits of this proposal, in terms of accuracy of the estimation and computational cost, are illustrated via numerical simulations.

Keywords:

state estimation; quantized data; Gaussian sum filtering; Gaussian sum smoothing; Gauss–Legendre quadrature

1. Introduction

It is well known that discrete-time dynamical systems can be described as first-order difference equations relating internal variables called states [1]. State estimation is a scientific discipline that studies methodologies and algorithms for estimating the state of dynamical systems from input–output measurements [2,3]. There are a variety of applications that use state estimation, such as control [4,5,6,7], parameter identification [8,9,10], power systems [11,12], fault detection [13,14,15,16,17], prognosis [18,19], cyber–physical systems [20], hydrologic and geophysical data assimilation [21,22], maritime tracking [23], consensus-based state estimation using wireless sensor networks [24,25,26], navigation systems [27], and transportation and highway traffic management [28,29,30], to mention a few. Depending on the measurements that are used, two algorithms of state estimation can be distinguished: filtering and smoothing. Filtering algorithms estimate the current state using measurements up to the current instant, and smoothing algorithms estimate the state at some time in the past using measurement up to the current instant [23,31].

In general, the experimental noisy data can suffer a significant loss of information introduced by low-resolution and cost-effective sensors [32] in the digitalization process [33]. Typically, the digitalization process encompasses a process known as quantization. Quantization is a nonlinear map that partitions the whole signal space and represents all of the values in each subspace by a single one [32]. In spite of the loss of information, the benefits of quantization have led to a number of applications in which quantized measurements arise. This occurs due to fundamental limitations on measuring equipment and bandwidth resources [34], digital and analog converters [35], and experimental designs where it is necessary to quantize the data in order to store it or minimize communication resource utilization [36]. In particular, estimation problems utilizing quantized measurements arise in networked control over limited-/finite-capacity communication channels, where usually, encoder–decoder state estimation schemes are used [37,38,39,40]. In addition, in order to deal with uncertain dynamic systems, robust estimation algorithms have also been developed; see, e.g., [37,38,41].

Currently, state estimation from quantized data has gained significant attention in a growing number of applications such as fault detection [42,43], networked control [42,44,45], and system identification [35,46,47,48]. For instance, in [49,50], Kalman-based state estimation algorithms were developed using multiple sensors for distributed systems. In [24], a Kalman smoothing algorithm was developed for any-time minimum-mean-squared error optimal-consensus-based state estimation using wireless sensor networks. In [27], an online smoothing algorithm was developed to estimate the positional and orientation parameters of integrated navigation systems utilizing a low-cost microelectromechanical system inertial sensor in near-real time.

In Figure 1, a usual representation of a process is shown, which is defined by the interconnection of three blocks: (1) an actuator, (2) a process, and (3) a sensor. This representation includes a link that can be a communication channel and a base station. The actuator input usually comes from a control system, and the output of the process is measured with noise by a sensor. The sensor introduces quantization to the noisy measurements. This representation has been used for state estimation and control in a microgrid with multiple distributed energy resources, where a dequantizer is used to reconstruct the received signal and then perform the standard Kalman filter [51]. In [42], a fault isolation filter for a discrete-time networked control system with multiple faults was developed using the Kalman filter, where the sensor measurements were transmitted only when the output signal was greater than a threshold. The authors in [52] dealt with a similar structure to the one in Figure 1 to estimate the vehicle sideslip angle of an in-vehicle networked system with sensor failure, dynamic quantization, and data dropouts. For more examples, see, e.g., [15,53].

For linear systems and Gaussian noises (without quantization), the optimal estimator of the system state is the celebrated Kalman filter [1] and smoother [22]. However, in the most general case, i.e., nonlinear systems and non-Gaussian noises, it is not possible to obtain an optimal estimator because the computation of some integrals in the filtering and smoothing equations is difficult or the integrals are just intractable. As mentioned above, the quantization process is a nonlinear map that results in a significant loss of information on the system dynamics, which produces a biased state estimation and incorrect characterization of the filtering and smoothing probability density functions (PDFs). In this context, several suboptimal filtering and smoothing algorithms for state-space systems with quantized data have been developed, for instance, standard- [49,54], unscented- [55], and extended- [36] Kalman filters for quantized data, in which some structural elements of the state-space models and the quantizer are exploited. Sequential Monte Carlo methods have been also used for filtering and smoothing with quantized data, where complex integrals are approximated by a set of weighted samples called particles [31], which define (approximately) a desired PDF. However, the most common difficulty in these methods is dealing with the quantizer model. Some approaches have been proposed for this purpose such as gradient-based approximation of the quantizer [56] or modeling the quantizer as uniformly distributed additive noise [32,57]. In [58], an approximation of the integral in the non-Gaussian probability mass function (PMF) of the quantized output given the state was proposed by using Gaussian quadrature rules in order to deal with binary data. This approximation naturally yields an explicit Gaussian mixture model (GMM) form for the PMF of the quantized output. In this paper, the approximation of this integral was extended by considering a more general (finite- and infinite-level) quantizer, and it was used to develop Gaussian sum filtering and smoothing algorithms.

Main Contributions

The main contributions of this work are:

Developing an explicit model (of the GMM form) for the PMF of the quantized output considering a finite- and infinite-level quantizer, to solve in closed-form filtering and smoothing recursions in a Bayesian framework;
Designing Gaussian sum filtering and smoothing algorithms to deal with quantized data, providing closed expressions for the state estimates and for the filtering and smoothing PDFs.

The filtering algorithm for quantized data presented in this paper includes, as a particular case, the filtering algorithm presented in [58], where only the case of binary data was considered. Additionally, the smoothing algorithm presented here, based on the approximation of the PMF of the quantized output given the state, is completely novel. The remainder of the paper is as follows: In Section 2, the problem of interest is defined. In Section 3, the characterization and approximation of the PMF of the quantized data given the state are presented, and using this explicit model, the Gaussian sum filtering and smoothing algorithms are developed. In Section 4, some examples are presented to show the benefits of this proposal in terms of accuracy and computational cost. Finally, in Section 5, conclusions are presented.

2. Statement of the Problem

2.1. System Model

Consider the state-space model for a discrete-time linear-time-invariant system with quantized output (see Figure 2):

\begin{matrix} (1) & x_{t + 1} & = A x_{t} + B u_{t} + w_{t}, \\ (2) & z_{t} & = C x_{t} + D u_{t} + v_{t}, \\ (3) & y_{t} & = q \{z_{t}\}, \end{matrix}

where

x_{t} \in R^{n}

is the state vector,

z_{t} \in R

is the nonquantized output,

y_{t} \in R

is the quantized output, and

u_{t} \in R^{m}

is the input of the system. The matrices are

A \in R^{n \times n}

B \in R^{n \times m}

C \in R^{1 \times n}

, and

D \in R^{1 \times m}

. The nonlinear map

q \{\cdot\}

is the quantizer. The state noise

w_{t} \in R^{n}

and the output noise

v_{t} \in R

are zero-mean white Gaussian noises with covariance matrix Q and R, respectively. The system in (1)–(2) can be described using the state transition and conditional nonquantized output probability distributions:

\begin{matrix} (4) & p (x_{1}) & = N_{x_{1}} (μ_{1}, P_{1}), \\ (5) & p (x_{t + 1} | x_{t}) & = N_{x_{t + 1}} (A x_{t} + B u_{t}, Q), \\ (6) & p (z_{t} | x_{t}) & = N_{z_{t}} (C x_{t} + D u_{t}, R), \end{matrix}

where

N_{x} (μ, P)

represents a PDF corresponding to a Gaussian distribution with mean

μ

and the covariance matrix P of the variable x. The initial condition

x_{1}

, the model noise

w_{t}

, and the measurement noise

v_{t}

are statistically independent random variables.

2.2. Quantizer Model

Let the quantizer

q \{\cdot\} : R \to V

be a map from the real line defined by the set of intervals

\{J_{i} \subset R : i \in I\}

to the finite or countable infinite set

V = \{β_{i} \in R : i \in I\}

[59], where

I

is a set of indices defined by the quantizer type. In this work, two quantizers were considered. The first an infinite-level quantizer, in which the output of the quantizer has infinite (countable) levels of quantization corresponding to the indices’ set:

\begin{matrix} I = \{\dots, 1, 2, \dots, L, \dots\} . \end{matrix}

(7)

The definition of the infinite-level quantizer is as follows (see Figure 3 (left)):

q \{z_{t}\} = \begin{matrix} β_{i} & if & z_{t} \in J_{i}, \end{matrix} i \in I,

(8)

where the sets

J_{i} = \{z_{t} : q_{i - 1} \leq z_{t} < q_{i}\}

are disjoint intervals and each

β_{i}

is the value that the quantizer takes in the region

J_{i}

. The second is a finite-level quantizer, in which the output of the quantizer is limited to a minimum and maximum values (saturated quantizer) corresponding to the indices’ set:

\begin{matrix} I = \{1, 2, \dots, L - 1, L\} . \end{matrix}

(9)

The definition of the finite-level quantizer (see Figure 3 (right)) is similar to (8), where

I

is defined in (9),

J_{1} = \{z_{t} : z_{t} < q_{1}\}

J_{L} = \{z_{t} : q_{L - 1} \leq z_{t}\}

, and

J_{i} = \{z_{t} : q_{i - 1} \leq z_{t} < q_{i}\}

with

i = 2 \dots, L - 1

2.3. Problem Definition

The problem of interest can be defined as follows: Given the available data

u_{1 : N} = \{u_{1}, u_{2}, \dots, u_{N}\}

and

y_{1 : N} = \{y_{1}, y_{2}, \dots, y_{N}\}

, where N is the data length, obtain the filtering and smoothing PDFs of the state given the quantized measurements,

p (x_{t} | y_{1 : t})

and

p (x_{t} | y_{1 : N})

, respectively the state estimators:

\begin{matrix} (10) & {\hat{x}}_{t | t} & = E \{x_{t} | y_{1 : t}\}, \\ (11) & {\hat{x}}_{t | N} & = E \{x_{t} | y_{1 : N}\}, \end{matrix}

and the corresponding covariance matrices of the estimation error:

\begin{matrix} (12) & Σ_{t | t} & = E \{(x_{t} - {\hat{x}}_{t | t}) {(x_{t} - {\hat{x}}_{t | t})}^{T} | y_{1 : t}\}, \\ (13) & Σ_{t | N} & = E \{(x_{t} - {\hat{x}}_{t | N}) {(x_{t} - {\hat{x}}_{t | N})}^{T} | y_{1 : N}\}, \end{matrix}

where

t \leq N

and

E \{x | y\}

denotes the conditional expectation of x given y.

3. Gaussian Sum Filtering and Smoothing for Quantized Data

Here, the Gaussian sum filtering and smoothing algorithms for quantized output data are explained in detail.

3.1. Gaussian Mixture Models

Gaussian mixture models refer to a convex combination of Gaussian densities corresponding to a random variable

ζ \in R^{n}

. Then, the PDF can be written as [60]:

p (ζ) = \sum_{i = 1}^{K} φ_{i} N_{ζ} (υ_{i}, Γ_{i})

subject to

φ_{i} > 0

and

\sum_{i = 1}^{K} φ_{i} = 1

, where

φ_{i}

is the ith mixing weight,

υ_{i} \in R^{n}

is the ith mean, and

Γ_{i} \in R^{n \times n}

is the ith covariance matrix. GMMs are used in a variety of applications to approximate non-Gaussian densities [61,62,63] and filtering and smoothing [64,65], to mention a few.

3.2. General Bayesian Framework

The well-known equations for Bayesian filtering (see, e.g., [31]) are given by:

\begin{matrix} (14) & p (x_{t} | y_{1 : t}) & = \frac{p (y_{t} | x_{t}) p (x_{t} | y_{1 : t - 1})}{p (y_{t} | y_{1 : t - 1})}, \\ (15) & p (x_{t + 1} | y_{1 : t}) & = \int p (x_{t + 1} | x_{t}) p (x_{t} | y_{1 : t}) d x_{t}, \end{matrix}

where

p (x_{t} | y_{1 : t})

and

p (x_{t + 1} | y_{1 : t})

are the measurement and time-update equations, respectively. The PDF

p (x_{t + 1} | x_{t})

is obtained from the model in (5), and

p (y_{t} | y_{1 : t - 1})

is a normalization constant. On the other hand, the well-known formula for Bayesian smoothing (see, e.g., [31]) is defined by:

p (x_{t} | y_{1 : N}) = p (x_{t} | y_{1 : t}) \int \frac{p (x_{t + 1} | y_{1 : N}) p (x_{t + 1} | x_{t})}{p (x_{t + 1} | y_{1 : t})} d x_{t + 1} .

(16)

However, the PDF in (16) is difficult to compute because of the division by a non-Gaussian density. This difficulty was overcome in [64] where the smoothing equations were separated in a two-fold formula filter: (i) the first formula, called the backward filter, defined by the following recursion:

\begin{matrix} (17) & p (y_{t + 1 : N} | x_{t}) & = \int p (y_{t + 1 : N} | x_{t + 1}) p (x_{t + 1} | x_{t}) d x_{t + 1}, \\ (18) & p (y_{t : N} | x_{t}) & = p (y_{t} | x_{t}) p (y_{t + 1 : N} | x_{t}), \end{matrix}

where

p (y_{t + 1 : N} | x_{t})

and

p (y_{t : N} | x_{t})

are the backward prediction and the backward- measurement-update equations, respectively, and (ii) the second formula, given by:

p (x_{t} | y_{1 : N}) = \frac{p (x_{t} | y_{1 : t - 1}) p (y_{t : N} | x_{t})}{p (y_{t : N} | y_{1 : t - 1})},

(19)

where

p (x_{t} | y_{1 : t - 1})

is the time-update equation from the filtering step,

p (y_{t : N} | y_{1 : t - 1})

is a normalization constant, and

p (y_{t : N} | x_{t})

is obtained using the backward recursion given in (17) and (18).

Due to the non-Gaussianity of

p (y_{t} | x_{t})

, the integrals in both the filtering and backward-filtering algorithms in (15) and (17), respectively, are difficult to compute or intractable. Widely used methods to deal with these integrals are Monte-Carlo-based algorithms, such as particle filtering and smoothing; see, e.g., [31,66]. These methods represent the posterior distributions

p (x_{t} | y_{1 : t})

and

p (x_{t} | y_{1 : N})

by a set of weighted random samples. However, in general, they exhibit a growing computational complexity as the model dimension increases. Here, an alternative method to compute the PMF

p (y_{t} | x_{t})

using Gauss–Legendre quadrature is proposed. This procedure results in a GMM form. Hence, the integrals in both the filtering and backward filtering in (15) and (17), respectively, can be computed in closed form.

Remark 1.

Notice that since

y_{t}

is a discrete random variable, the measurement-update equation in (14) and the backward-measurement-update equation in (18) comprise PDFs and a PMF. Hence, generalized probability density functions are used here; see, e.g., [67].

3.3. Computing an Explicit Model of $p (y_{t} | x_{t})$

From (8), it is observed that the output

z_{t} \in J_{i}

is mapped to a single value

β_{i}

. Then, the probability that

y_{t}

takes the value

β_{i}

is the same as the probability of

z_{t}

belonging to set

J_{i}

, as shown in Figure 4. Notice that the quantizer regions

J_{i}

include finite and semi-infinite intervals.

In the following theorem, the characterization of

p (y_{t} | x_{t})

is formalized via the computation of the probability

P \{z_{t} \in J_{i}\}

shown in Figure 4 through the integral definition of probabilities [68]. Therefore, the approximation of this integral by using the Gauss–Legendre quadrature rule is presented (see, e.g., [69]).

Theorem 1.

Consider the system (1)–(3), the infinite- and finite-level quantizers defined in (8). Then, the PMF of the discrete random variable

y_{t}

given the state

x_{t}

is given by:

p (y_{t} | x_{t}) = \int_{a_{t}}^{b_{t}} N_{v_{t}} (0, R) d v_{t},

(20)

where

a_{t}

and

b_{t}

are functions of the boundary values of each region of the quantizers and are defined in Table 1. In addition, the integral in (20) can be approximated using the Gauss–Legendre quadrature rule, yielding:

p (y_{t} | x_{t}) \approx \sum_{τ = 1}^{K} ς_{t}^{τ} N_{η_{t}^{τ}} (C x_{t} + D u_{t} + μ_{t}^{τ}, R),

(21)

where K is the number of points from the Gauss–Legendre quadrature rule (defined by the user),

ς_{t}^{τ}

η_{t}^{τ}

, and

μ_{t}^{τ}

are defined in Table 1, and

ω_{τ}

and

ψ_{τ}

are weights and points defined by the quadrature rule, given in, e.g., [69].

Proof.

From the infinite- and finite-level quantizers, it is observed that the random variable

y_{t}

can only take the values

β_{i}

with i in the indices’ sets given in (7) and (9). Then, the probability of

y_{t} = \dots, β_{1}, β_{2}, \dots, β_{L}, \dots

y_{t} = β_{1}, β_{2}, \dots, β_{L - 1}, β_{L}

is the same as the probability that the random variable

z_{t}

belongs to the sets

J_{i}

. This probability can be obtained from the distribution function as follows:

P \{y_{t} = β_{i} | x_{t}\} = P \{z_{t} \in J_{i} | x_{t}\},

(22)

where

P \{\cdot\}

denotes probability. Considering the infinite-level quantizer and the output equation in (2), the following expressions are obtained for

y_{t} = β_{i}

with

i = \dots, 1, 2 \dots, L, \dots

\begin{matrix} P \{y_{t} = β_{i} | x_{t}\} & = P \{q_{i - 1} \leq z_{t} < q_{i} | x_{t}\}, \\ = P \{a_{t} \leq v_{t} < b_{t} | x_{t}\}, \end{matrix}

(23)

where

a_{t} = q_{i - 1} - C x_{t} - D u_{t}

and

b_{t} = q_{i} - C x_{t} - D u_{t}

. Additionally, for the finite-level quantizer, (23) holds for

y_{t} = β_{i}

with

i = 2, \dots, L - 1

, and for

y_{t} = β_{1}

and

y_{t} = β_{L}

, the following holds:

\begin{matrix} (24) & P \{v_{t} < b_{t} | x_{t}\} & if y_{t} = β_{1}, \\ (25) & P \{v_{t} \geq a_{t} | x_{t}\} & if y_{t} = β_{L}, \end{matrix}

where

b_{t} = q_{1} - C x_{t} - D u_{t}

and

a_{t} = q_{L - 1} - C x_{t} - D u_{t}

. Then, using the fact that

p (v_{t}) = N_{v_{t}} (0, R)

and using (23)–(25), the integral in (20) can be obtained, where the integral limits are given in Table 1.

On the other hand, the Gauss–Legendre quadrature for approximating a definite integral over the interval

[- 1, 1]

is given by:

\int_{- 1}^{1} f (r) d r \approx \sum_{τ = 1}^{K} ω_{τ} f (ψ_{τ}),

(26)

where

ω_{τ}

and

ψ_{τ}

are the quadrature weights and the roots of the order K Legendre polynomial, respectively; see, e.g., [69]. Notice that the integral over

[a_{t}, b_{t}]

can be mapped onto the interval

[- 1, 1]

using:

\int_{a_{t}}^{b_{t}} f (r) d r = \frac{b_{t} - a_{t}}{2} \int_{- 1}^{1} f (\frac{b_{t} - a_{t}}{2} r + \frac{b_{t} + a_{t}}{2}) d r,

(27)

and using the definition in (26), this integral is approximated by:

\int_{a_{t}}^{b_{t}} f (r) d r \approx \frac{b_{t} - a_{t}}{2} \sum_{τ = 1}^{K} ω_{τ} f (\frac{b_{t} - a_{t}}{2} ψ_{τ} + \frac{b_{t} + a_{t}}{2}) .

(28)

Defining

a_{t} = q_{i - 1} - C x_{t} - D u_{t}

b_{t} = q_{i} - C x_{t} - D u_{t}

, and

f (v_{t}) = N_{v_{t}} (0, R)

with

i = \dots, 1, 2 \dots, L, \dots

, the approximation of

p (y_{t} | x_{t})

given in (21) for the infinite-level quantizer is derived.

For the finite-level quantizer, (28) holds for

y_{t} = β_{i}

with

i = 2, \dots, L - 1

. Then, the approximation of

p (y_{t} | x_{t})

for

y_{t} = β_{1}

and

y_{t} = β_{L}

is defined. First, it is observed that the integral over the semi-infinite interval

[a_{t}, \infty)

can be mapped onto the interval

(0, 1]

using

r = a_{t} + (1 - ξ) / ξ

, so that:

\int_{a_{t}}^{\infty} f (r) d r = \int_{0}^{1} f (a_{t} + \frac{1 - ξ}{ξ}) \frac{d ξ}{ξ^{2}} .

(29)

Then, using an appropriate change of variable, it can be mapped onto the interval

(- 1, 1]

, yielding:

\int_{a_{t}}^{\infty} f (r) d r = \int_{- 1}^{1} f (a_{t} + \frac{1 - s}{1 + s}) \frac{2 d s}{{(1 + s)}^{2}} .

(30)

Using the Gauss–Legendre approximation given in (26), the approximation of

p (y_{t} | x_{t})

for

y_{t} = β_{L}

is obtained as follows:

\int_{a_{t}}^{\infty} f (r) d r \approx \sum_{τ = 1}^{K} ω_{τ} f (a_{t} + \frac{1 - ψ_{τ}}{1 + ψ_{τ}}) \frac{2}{{(1 + ψ_{τ})}^{2}} .

(31)

Defining

a_{t} = q_{L - 1} - C x_{t} - D u_{t}

b = \infty

, and

f (v_{t}) = N_{v_{t}} (0, R)

, (21) is obtained. A similar procedure for the integral over the semi-infinite interval

(- \infty, b_{t}]

can be applied using

r = b_{t} - (1 - ξ) / ξ

to find the approximation of

p (y_{t} | x_{t})

for

y_{t} = β_{1}

, as follows:

\int_{- \infty}^{b_{t}} f (r) d r \approx \sum_{τ = 1}^{K} ω_{τ} f (b_{t} - \frac{1 - ψ_{τ}}{1 + ψ_{τ}}) \frac{2}{{(1 + ψ_{τ})}^{2}} .

(32)

Defining

a_{t} = - \infty

b_{t} = q_{1} - C x_{t} - D u_{t}

, and

f (v_{t}) = N_{v_{t}} (0, R)

, (21) is obtained. This completes the proof. □

Remark 2.

Notice that any quadrature rule, such as Newton–Cotes, Gauss–Laguerre, and Gauss–Hermite (see, e.g., [69,70]), used to approximate the integral in (20), yields a weighted sum of Gaussian distributions evaluated as a linear function of the values

q_{i}

and the state

x_{t}

. Furthermore, it is possible to interpret

p (y_{t} | x_{t})

as a weighted sum of Gaussian distributions in the random variable

x_{t}

. This weighted sum is denoted as the GMM structure. Thus, this structure is considered for developing the Gaussian sum filtering and smoothing algorithms that deal with quantized data.

Remark 3.

Notice that in [70], a suboptimal filtering algorithm called the quadrature Kalman filter, where a linear regression is used to linearize the nonlinear process and measurement functions using the Gauss–Hermite quadrature rule, was developed. This approach is not directly applicable to the problem of interest in this paper, so that the quantizer is a nondifferentiable nonlinearity.

On the other hand, the cubature Kalman filter [70] and smoother [71] are approaches that use the spherical–radial cubature rule to approximate the n-dimensional integrals when computing the expected values in (10)–(13) in a nonlinear state-space model. These integral approximations are obtained under the assumption that

p (x_{t} | y_{1 : t})

and

p (x_{t} | y_{1 : N})

are Gaussian distributed. The difference between these approaches and the proposed method in this paper is that the Gauss quadrature rule was used to approximate the integral in the probabilistic model of

p (y_{t} | x_{t})

. It is clear that in this paper, the Gaussian assumption in the filtering and smoothing densities was not used. In fact, it is shown that

p (x_{t} | y_{1 : t})

and

p (x_{t} | y_{1 : N})

are non-Gaussian PDFs.

3.4. Understanding the Gaussian Sum Filtering and Smoothing Algorithms

The general Bayesian framework for filtering leads to the product

p (y_{t} | x_{t})

and

p (x_{t} | y_{1 : t - 1})

, as shown in (14). From the definition of

p (y_{t} | x_{t})

in (21), it is observed that

p (x_{t} | y_{1 : t})

results in the product of two GMMs. This, in turn, results in a GMM with more components than the individual factors

p (y_{t} | x_{t})

and

p (x_{t} | y_{1 : t - 1})

. This implies that the time-update equation

p (x_{t + 1} | y_{1 : t})

in (15) is also a GMM distribution. A similar situation occurs in the backward-measurement-update equation

p (y_{t : N} | x_{t})

in (18), which is the product of two GMM structures

p (y_{t} | x_{t})

and

p (y_{t + 1 : N} | x_{t})

. This yields another GMM structure with more components than the individual factors, which implies that the backward-prediction equation in (17) is also a GMM structure. For the clarity of the presentation, reducing the product of two GMMs (structures) into one is necessary.

In order to understand the transformation of the product of two summations into another summation, the following sums are considered:

g = \sum_{τ = 1}^{K} g_{τ}

and

f = \sum_{ℓ = 1}^{M} f_{ℓ}

. Then, for each two-tuple

(τ, ℓ)

where

τ = 1, \dots, K

and

ℓ = 1, \dots, M

, the product

h = f g

is another summation and has

K M

terms indexed by

k = (ℓ - 1) K + τ

. Then, reordering these terms, the following sum is obtained:

h = \sum_{k = 1}^{K M} h_{k}

, where

h_{k} = f_{ℓ} g_{τ}

3.5. Gaussian Sum Filtering for Quantized Data

Using the approximation of the PMF

p (y_{t} | x_{t})

defined in Theorem 1, the Gaussian sum filtering algorithm can be derived using (14) and (15) as follows:

Theorem 2 (Gaussian sum filter).

Consider the system in (4)–(6) and the approximation of

p (y_{t} | x_{t})

in (21). The filtering equations for state-space systems with quantized output data are defined as follows:

Initialization: For

t = 1

, the predicted distribution is

p (x_{1}) = N_{x_{1}} (μ_{1}, P_{1})

, where

p (x_{1})

is the prior distribution of the initial state. Then, for

t = 1, \dots, N

, the measurement-update and the time-update equations for the Gaussian sum filtering are defined as follows:

Measurement update: The PDF of the current state

x_{t}

given the current and previous measurements, that is

y_{1}, \dots, y_{t}

, is the GMM given by:

p (x_{t} | y_{1 : t}) = \sum_{k = 1}^{M_{t | t}} γ_{t | t}^{k} N_{x_{t}} ({\hat{x}}_{t | t}^{k}, Σ_{t | t}^{k}),

(33)

where

γ_{t | t}^{k}

{\hat{x}}_{t | t}^{k}

, and

Σ_{t | t}^{k}

are given in Appendix B,

M_{t | t} = K M_{t | t - 1}

, and

M_{t | t - 1}

is the number of Gaussian components in the time update step. The initial values satisfy that

M_{1 | 0} = 1

γ_{1 | 0} = 1

{\hat{x}}_{1 | 0} = μ_{1}

, and

Σ_{1 | 0} = P_{1}

Time update: The PDF of the future state

x_{t + 1}

, one-step-ahead prediction given the measurements

y_{1}, \dots, y_{t}

, is the GMM given by:

p (x_{t + 1} | y_{1 : t}) = \sum_{k = 1}^{M_{t + 1 | t}} γ_{t + 1 | t}^{k} N_{x_{t + 1}} ({\hat{x}}_{t + 1 | t}^{k}, Σ_{t + 1 | t}^{k}),

(34)

where

M_{t + 1 | t} = M_{t | t}

γ_{t + 1 | t}^{k}

{\hat{x}}_{t + 1 | t}^{k}

, and

Σ_{t + 1 | t}^{k}

are given in Appendix B.

Proof.

Consider the recursion in (14) and (15). The required PDF

p (x_{t + 1} | x_{t})

and PMF

p (y_{t} | x_{t})

are obtained from (5) and (21), respectively. Then, using the measurement-update equations in (14) and Lemma 2 in [65] (p. 86), the following expression is obtained:

p (x_{t} | y_{1 : t}) \propto \sum_{τ = 1}^{K} \sum_{ℓ = 1}^{M_{t | t - 1}} ς_{t}^{τ} γ_{t | t - 1}^{ℓ} N_{η_{t}^{τ}} (κ_{t}^{ℓ τ}, V_{t}^{ℓ}) N_{x_{t}} ({\hat{x}}_{t | t}^{k}, Σ_{t | t}^{k})

(35)

where

κ_{t}^{ℓ τ}

V_{t}^{ℓ}

{\hat{x}}_{t | t}^{k}

, and

Σ_{t | t}^{k}

are defined in (A11), (A12), (A7), and (A8), respectively. Notice that

N_{η_{t}^{τ}} (κ_{t}^{ℓ τ}, V_{t}^{ℓ})

is a coefficient since the measurement

y_{t}

is available. Then, rewriting the double summation in (35) as a single one with a new index

k = (ℓ - 1) K + τ

, the measurement-update equation in (33) is derived defining

M_{t | t} = K M_{t | t - 1}

, the weights

{\bar{γ}}_{t | t}^{k}

as in (A9), and computing the corresponding normalization constant. On the other hand, using the time-update equation in (15) and Lemma 2 in [65] (p. 86), the following equation is obtained:

p (x_{t + 1} | y_{1 : t}) = \sum_{k = 1}^{M_{t | t}} γ_{t | t}^{k} \int N_{x_{t}} (m_{t}^{k}, S_{t}^{k}) N_{x_{t + 1}} ({\hat{x}}_{t + 1 | t}^{k}, Σ_{t + 1 | t}^{k}) d x_{t},

(36)

where

{\hat{x}}_{t + 1 | t}^{k}

and

Σ_{t + 1 | t}^{k}

are defined in (A14), and (A15), respectively. Then, solving this integral, the time-update equation in (34) can be derived defining

M_{t + 1 | t} = M_{t | t}

, the weights

γ_{t + 1 | t}^{k}

as in (A13), and computing the corresponding normalization constant. This completes the proof. □

3.6. Computing the State Estimator ${\hat{x}}_{t | t}$ from a GMM

Provided

p (x_{t} | y_{1 : t})

in (33) as a GMM, the state estimator given in (10) and the covariance matrix of the estimation error in (12) can be computed as follows:

\begin{matrix} (37) & {\hat{x}}_{t | t} & = \sum_{k = 1}^{M_{t | t}} γ_{t | t}^{k} {\hat{x}}_{t | t}^{k}, \\ (38) & Σ_{t | t} & = \sum_{k = 1}^{M_{t | t}} γ_{t | t}^{k} [Σ_{t | t}^{k} + ({\hat{x}}_{t | t}^{k} - {\hat{x}}_{t | t}) {({\hat{x}}_{t | t}^{k} - {\hat{x}}_{t | t})}^{T}] . \end{matrix}

3.7. Backward Filtering for Quantized Data

Using the approximation of the PMF

p (y_{t} | x_{t})

defined in Theorem 1, the backward filter recursion can be derived as follows:

Theorem 3 (Backward filtering).

Consider the system in (4)–(6) and the approximation of

p (y_{t} | x_{t})

in (21). Then, the backward filter for state-space systems with quantized output data is:

Initialization: For

t = N

, the backward measurement update is given by:

\begin{matrix} p (y_{N} | x_{N}) = \sum_{k = 1}^{S_{N | N}} ϵ_{N | N}^{k} λ_{N | N}^{k} exp \{- \frac{1}{2} (x_{N}^{T} F_{N | N}^{k} x_{N} - 2 G_{N | N}^{k T} x_{N} + H_{N | N}^{k})\}, \end{matrix}

(39)

where

S_{N | N} = K

, and:

\begin{matrix} ϵ_{N | N}^{k} & = ς_{N}^{k}, \\ F_{N | N}^{k} & = C^{T} R^{- 1} C, \end{matrix} \begin{matrix} λ_{N | N}^{k} & = {(det \{2 π R\})}^{- 1 / 2} \\ G_{N | N}^{k T} & = θ_{N}^{k T} R^{- 1} C, \end{matrix} \begin{matrix} θ_{N}^{k} & = η_{N}^{k} - D u_{N} - μ_{N}^{k}, \\ H_{N | N}^{k} & = θ_{N}^{k T} R^{- 1} θ_{N}^{k}, \end{matrix}

(40)

where

ς_{N}^{k}

η_{N}^{k}

, and

μ_{N}^{k}

are defined in Table 1. Then, for

t = N - 1, \dots, 1

, the backward prediction and the backward-measurement-update equations are defined as follows:

Backward predictions: The backward-prediction equation is given by:

\begin{matrix} p (y_{t + 1 : N} | x_{t}) = \sum_{k = 1}^{S_{t | t + 1}} ϵ_{t | t + 1}^{k} λ_{t | t + 1}^{k} exp \{- \frac{1}{2} (x_{t}^{T} F_{t | t + 1}^{k} x_{t} - 2 G_{t | t + 1}^{k T} x_{t} + H_{t | t + 1}^{k})\}, \end{matrix}

(41)

where

S_{t | t + 1} = S_{t + 1 | t + 1}

S_{t + 1 | t + 1}

is the number of components in the backward-measurement-update step, and

ϵ_{t | t + 1}^{k}

λ_{t | t + 1}^{k}

F_{t | t + 1}^{k}

G_{t | t + 1}^{k T}

H_{t | t + 1}^{k}

are given in Appendix C.

Backward-measurement update: The distribution of

y_{t : N} | x_{t}

evaluated at

y_{t}, \dots, y_{N}

is:

\begin{matrix} p (y_{t : N} | x_{t}) = \sum_{k = 1}^{S_{t | t}} ϵ_{t | t}^{k} λ_{t | t}^{k} exp \{- \frac{1}{2} (x_{t}^{T} F_{t | t}^{k} x_{t} - 2 G_{t | t}^{k T} x_{t} + H_{t | t}^{k})\}, \end{matrix}

(42)

where

S_{t | t} = K S_{t | t + 1}

S_{t | t + 1}

is the number of components in the backward-prediction step, and

ϵ_{t | t}^{k}

λ_{t | t}^{k}

F_{t | t}^{k}

G_{t | t}^{k T}

H_{t | t}^{k}

are given in Appendix C.

Proof.

Consider the recursion in (17) and (18). The required PDF

p (x_{t + 1} | x_{t})

and PMF

p (y_{t} | x_{t})

are given by Equations (5) and (21), respectively. The proof is carried out by induction in reverse time. First, it is verified that it holds for

t = N - 1

, then it is assumed to be true for

t = s + 1

, and finally, it is verified that it holds for

t = s

. Notice that the recursion starts in

t = N

with

p (y_{N : N} | x_{N}) = p (y_{N} | x_{N})

, which is defined in (39). From (17), at time

t = N - 1

, the backward-prediction step is defined as:

p (y_{N : N} | x_{N - 1}) = \int p (y_{N : N} | x_{N}) p (x_{N} | x_{N - 1}) d x_{N},

(43)

where

p (y_{N : N} | x_{N})

is given by (39) and

p (x_{N} | x_{N - 1})

is defined by the system model in (5). From the definition given in (43), the equation below is obtained:

\begin{matrix} p (y_{N : N} | x_{N - 1}) = \sum_{k = 1}^{S_{N | N}} \frac{ϵ_{N | N}^{k} λ_{N | N}^{k}}{\sqrt{det \{2 π Q\}}} \int exp \{- \frac{1}{2} [x_{N}^{T} {\bar{P}}_{N}^{k} x_{N} - 2 {\bar{V}}_{N}^{k T} x_{N} + {\bar{S}}_{N}^{k}]\} d x_{N}, \end{matrix}

(44)

where

{\bar{P}}_{N}^{k} = (F_{N | N}^{k} + Q^{- 1})

{\bar{V}}_{N}^{k} = (G_{N | N}^{k} + J_{N - 1})

{\bar{S}}_{N}^{k} = (H_{N | N}^{k} + L_{N - 1})

J_{N - 1}^{T} = {(A x_{N - 1} + B u_{N - 1})}^{T} Q^{- 1}

, and

L_{N - 1} = {(A x_{N - 1} + B u_{N - 1})}^{T} Q^{- 1} (A x_{N - 1} + B u_{N - 1})

. Then, completing the square and solving the integral in (44), it is obtained:

\begin{matrix} p (y_{N : N} | x_{N - 1}) & = \sum_{k = 1}^{S_{N - 1 | N}} ϵ_{N - 1 | N}^{k} λ_{N - 1 | N}^{k} \\ \times exp \{- \frac{1}{2} (x_{N - 1}^{T} F_{N - 1 | N}^{k} x_{N - 1} - 2 G_{N - 1 | N}^{k T} x_{N - 1} + H_{N - 1 | N}^{k})\}, \end{matrix}

(45)

where

S_{N - 1 | N} = S_{N | N}

and the remaining terms in (45) are defined in (A22)–(A26), but evaluated at

t = N - 1

. The backward-measurement-update step in (18) is as follows:

p (y_{N - 1 : N} | x_{N - 1}) = p (y_{N - 1} | x_{N - 1}) p (y_{N : N} | x_{N - 1}) .

(46)

Thus, using (45) in the definition given in (46), it is obtained:

\begin{matrix} p (y_{N - 1 : N} | x_{N - 1}) & = \sum_{τ = 1}^{K} \sum_{ℓ = 1}^{S_{N - 1 | N}} \frac{ς_{N - 1}^{τ} ϵ_{N - 1 | N}^{ℓ} λ_{N - 1 | N}^{ℓ}}{\sqrt{det \{2 π R\}}} \\ \times exp \{- \frac{1}{2} (x_{N - 1}^{T} {\tilde{F}}_{N - 1} x_{N - 1} - 2 {\tilde{G}}_{N - 1}^{τ T} x_{N - 1} + {\tilde{H}}_{N - 1}^{τ})\} \\ \times exp \{- \frac{1}{2} (x_{N - 1}^{T} F_{N - 1 | N}^{ℓ} x_{N - 1} - 2 G_{N - 1 | N}^{ℓ T} x_{N - 1} + H_{N - 1 | N}^{ℓ})\}, \end{matrix}

(47)

where

{\tilde{F}}_{N - 1} = C^{T} R^{- 1} C

{\tilde{G}}_{N - 1}^{τ T} = θ_{N - 1}^{τ T} R^{- 1} C

, and

{\tilde{H}}_{N - 1}^{τ} = θ_{N - 1}^{τ T} R^{- 1} θ_{N - 1}^{τ}

with

θ_{N - 1}^{τ} = η_{N - 1}^{τ} - D u_{N - 1} - μ_{N - 1}^{τ}

. Finally, rewriting the double summation in the last equation into a single one with a new index

k = (ℓ - 1) K + τ

results in:

\begin{matrix} p (y_{N - 1 : N} | x_{N - 1}) = \sum_{k = 1}^{S_{N - 1 | N - 1}} ϵ_{N - 1 | N - 1}^{k} λ_{N - 1 | N - 1}^{k} \\ \times exp \{- \frac{1}{2} [x_{N - 1}^{T} F_{N - 1 | N - 1}^{k} x_{N - 1} - 2 G_{N - 1 | N - 1}^{k T} x_{N - 1} + H_{N - 1 | N - 1}^{k}]\}, \end{matrix}

(48)

where

S_{N - 1 | N - 1} = K S_{N - 1 | N}

and the remaining terms in (48) are the same as the ones shown in (A16)–(A21), but evaluated at

t = N - 1

. Thus, it is concluded that it holds for

t = N - 1

. A similar procedure was applied to find that, for both the backward-prediction and backward-measurement-update steps in the backward filter, it yields the same expressions in (A22)–(A26) and (A16)–(A21), but evaluated at

t = s

. Thus, it is concluded that Theorem 3 holds for all t.

From Theorems 2 and 3, it is clear that the number of elements in the mixture grows exponentially with every iteration, making the algorithm computationally intractable after a few iterations. In addition, it would be necessary to save and manage a large amount of information in every iteration of the Gaussian sum filtering and in the backward recursion. Therefore, an algorithm that reduces the number of GMM components should be implemented in every iteration of these two algorithms in order to keep the number of components bounded. Different methods have been proposed to perform this kind of procedure, termed Gaussian sum reduction, such us pruning, joining, and integral-squared-error-based methods; see [70]. In this work, the Kullback–Leibler approach for Gaussian sum reduction proposed in [64,72] was used. The idea behind the Gaussian sum reduction is to transform the GMM

{\{φ_{i}, υ_{i}, Γ_{i}\}}_{i = 1}^{J}

into a GMM

{\{φ_{i}, υ_{i}, Γ_{i}\}}_{i = 1}^{S}

, where

1 \leq S \leq J

. In [64], it was suggested to use a measure of dissimilarity between two components and pooling the pair of components that minimize this measure. Then, based on this idea, in [72], the Kullback–Leibler information number was used as the measure of dissimilarity. The author in [72] provided an algorithm to merge two components so that the merged component preserves the first- and second-order moments of the original two components, which is given by:

(φ_{i j}, υ_{i j}, Γ_{i j}) = M \{(φ_{i}, υ_{i}, Γ_{i}), (φ_{j}, υ_{j}, Γ_{j})\},

(49)

where

M \{\cdot, \cdot\}

is a merging function such that:

\begin{matrix} (50) & φ_{i j} & = φ_{i} + φ_{j}, \\ (51) & υ_{i j} & = {\bar{φ}}_{i} υ_{i} + {\bar{φ}}_{j} υ_{j}, \\ (52) & Γ_{i j} & = {\bar{φ}}_{i} Γ_{i} + {\bar{φ}}_{j} Γ_{j} + {\bar{φ}}_{i} {\bar{φ}}_{j} (υ_{i} - υ_{j}) {(υ_{i} - υ_{j})}^{T}, \end{matrix}

where

{\bar{φ}}_{i} = φ_{i} / (φ_{i} + φ_{j})

and

{\bar{φ}}_{j} = φ_{j} / (φ_{i} + φ_{j})

. The merging function applied to two components i and j (

i \neq j

) minimizes the dissimilarity measure

D (i, j)

, defined as:

\begin{matrix} D (i, j) = \frac{1}{2} [φ_{i j} log det \{Γ_{i j}\} - φ_{i} log det \{Γ_{i}\} - φ_{j} log det \{Γ_{j}\}], \end{matrix}

(53)

The function

D (i, j)

satisfies

D (i, j) = D (j, i)

and

D (i, i) = 0

. This implies that the total number of combinations to merge is reduced to

0.5 J (J - 1)

. The authors in [73] used Runnalls’ algorithm in a Bayesian filtering environment, and they modified it to include a user-defined threshold for the number of components after reduction and a user-defined threshold

ε

that satisfies

D (i, j) < ε

. In Algorithm 1, the steps for implementing the Gaussian sum filtering are summarized.

Algorithm 1 Gaussian sum filter algorithm for quantized output data

₁ Input: The PDF of the initial state

p (x_{1})

, e.g.,

M_{1 | 0} = 1

γ_{1 | 0} = 1

{\hat{x}}_{1 | 0} = μ_{1}

Σ_{1 | 0} = P_{1}

. The points of the Gauss–Legendre quadrature

{\{ω_{τ}, ψ_{τ}\}}_{τ = 1}^{K}

₂₀ Output: The state estimation in (10), the covariance matrix of the estimation error

in (12), the filtering PDFs

p (x_{t} | y_{1 : t})

, the predictive PDFs

p (x_{t + 1} | y_{1 : t})

, and the set

\{ς_{t}^{τ}, η_{t}^{τ}, μ_{t}^{τ}\}

, for

t = 1, \dots, N

The backward recursion in Theorem 3 is used to obtain the smoothing PDF in (19). For this purpose,

p (y_{t : N} | x_{t})

is converted into a GMM structure of the random variable

x_{t}

. Then, the Gaussian sum reduction algorithm is applied to the GMM structure of

p (y_{t : N} | x_{t})

to obtain:

p (y_{t : N} | x_{t}) = \sum_{k = 1}^{S_{red}} δ_{t | t}^{k} N_{x_{t}} (z_{t | t}^{k}, U_{t | t}^{k}),

(54)

where

S_{red}

is the number of Gaussian components kept after the Gaussian reduction procedure and

δ_{t | t}^{k}

z_{t | t}^{k}

, and

U_{t | t}^{k}

are the corresponding weight, mean, and covariance matrix. This reduced GMM structure is used to obtain the smoothing PDFs. However, for the next recursion in the backward filter, reconverting the reduced GMM structure into the backward filter form to obtain the reduced version of the backward-measurement-update step in (42) is required. This conversion process between the backward filter and GMM structure is summarized in Lemma A3 in Appendix A. In Algorithm 2, the steps for implementing the backward filter are summarized.

Algorithm 2 Backward-filtering algorithm for quantized output data

₁ Input: The initial backward measurement

p (y_{N} | x_{N})

given in (39), e.g.,

S_{N | N}

ϵ_{N | N}^{k}

λ_{N | N}^{k}

F_{N | N}^{k}

G_{N | N}^{k T}

, and

H_{N | N}^{k}

. The set

\{ς_{t}^{τ}, η_{t}^{τ}, μ_{t}^{τ}\}

for

t = 1, \dots, N

computed in Algorithm 1.

₁₉ Output: The backward prediction

p (y_{t + 1 : N} | x_{t})

and the backward measurement

update

p (y_{t : N} | x_{t})

for

t = N, \dots, 1

3.8. Smoothing Algorithm with Quantized Data

In order to obtain the smoothing PDF in (19), the GMM structure

p (y_{t : N} | x_{t})

in the random variable

x_{t}

given in (54) is used. This GMM structure is multiplied by the time-update equation

p (x_{t} | y_{1 : t - 1})

of the filtering algorithm given in (34). Then, the smoothing PDF is obtained from the following:

Theorem 4 (Gaussian sum smoothing).

Consider the system in (4)–(6). Given

p (x_{t} | y_{1 : t - 1})

p (x_{N} | y_{1 : N})

and

p (y_{t : N} | x_{t})

, then the smoothing PDF at time

t = N

is given by

p (x_{N} | y_{1 : N})

, and for

t = N - 1, \dots, 1

, the PDF

p (x_{t} | y_{1 : N})

is a GMM given by:

p (x_{t} | y_{1 : N}) = \sum_{k = 1}^{S_{t | N}} ϵ_{t | N}^{k} N_{x_{t}} ({\hat{x}}_{t | N}^{k}, Σ_{t | N}^{k}),

(55)

where

S_{t | N} = M_{t | t - 1} S_{r e d}

M_{t | t - 1}

and

S_{r e d}

are given in (34) and (54), respectively, and

ϵ_{t | N}^{k}

{\hat{x}}_{t | N}^{k}

, and

Σ_{t | N}^{k}

are given in Appendix D.

Proof.

Consider the definition of Bayesian smoothing given in (19), the time-update equation

p (x_{t} | y_{1 : t - 1})

obtained from (34), and the reduced version of the measurement-update step in the backward filter

p (y_{t : N} | x_{t})

defined in (54). For

t = N - 1, N - 2, \dots, 1

, it is obtained:

p (x_{t} | y_{1 : N}) \propto \sum_{τ = 1}^{M_{t | t - 1}} \sum_{ℓ = 1}^{S_{red}} γ_{t | t - 1}^{τ} δ_{t | t}^{ℓ} N_{x_{t}} (z_{t | t}^{ℓ}, U_{t | t}^{ℓ}) N_{x_{t}} ({\hat{x}}_{t | t - 1}^{τ}, Σ_{t | t - 1}^{τ}) .

(56)

Defining

g (x_{t}) = N_{x_{t}} (z_{t | t}^{ℓ}, U_{t | t}^{ℓ}) N_{x_{t}} ({\hat{x}}_{t | t - 1}^{τ}, Σ_{t | t - 1}^{τ})

, the following equation is derived:

\begin{matrix} g (x_{t}) = \frac{exp \{- \frac{1}{2} [ϕ_{1 t}^{ℓ} + ϕ_{2 t}^{τ} - ϕ_{3 t}^{ℓ τ}]\} N_{x_{t}} ({(L_{t}^{ℓ τ})}^{- 1} ρ_{t}^{ℓ τ}, {(L_{t}^{ℓ τ})}^{- 1})}{{(2 π)}^{\frac{n}{2}} \sqrt{det \{L_{t}^{ℓ τ}\} det \{U_{t | t}^{ℓ}\} det \{Σ_{t | t - 1}^{τ}\}}}, \end{matrix}

(57)

where

L_{t}^{ℓ τ}

ρ_{t}^{ℓ τ}

ϕ_{1 t}^{ℓ}

ϕ_{2 t}^{τ}

, and

ϕ_{3 t}^{ℓ τ}

are defined in (A32), (A31), (A34), (A35), and (A33), respectively. Next, expressing the double summation in (56) as a single one with a new index

k = (ℓ - 1) M_{t | t - 1} + τ

, it is obtained:

\begin{matrix} p (x_{t} | y_{1 : N}) \propto \sum_{k = 1}^{S_{t | N}} {\bar{ϵ}}_{t | N}^{k} N_{x_{t}} ({\hat{x}}_{t | N}^{k}, Σ_{t | N}^{k}), \end{matrix}

(58)

where

S_{t | N} = M_{t | t - 1} S_{red}

and

{\bar{ϵ}}_{t | N}^{k}

Σ_{t | N}^{k}

, and

{\hat{x}}_{t | N}^{k}

are given in (A30), (A28), and (A29), respectively. Finally, the smoothing PDF in (55) is obtained by computing the normalized wights as (A27), and this completes the proof. □

Provided

p (x_{t} | y_{1 : N})

in (55) as a GMM, the state estimator given in (11) and the covariance matrix of the estimation error in (13) can be computed as follows:

\begin{matrix} (59) & {\hat{x}}_{t | N} & = \sum_{k = 1}^{S_{t | N}} ϵ_{t | N}^{k} {\hat{x}}_{t | N}^{k}, \\ (60) & Σ_{t | N} & = \sum_{k = 1}^{S_{t | N}} ϵ_{t | N}^{k} [Σ_{t | N}^{k} + ({\hat{x}}_{t | N}^{k} - {\hat{x}}_{t | N}) {({\hat{x}}_{t | N}^{k} - {\hat{x}}_{t | N})}^{T}] . \end{matrix}

In Algorithm 3, the steps to implement the Gaussian sum smoothing are summarized.

Algorithm 3 Gaussian sum smoothing algorithm for quantized output data

₁ Input: The PDFs

p (x_{t} | y_{1 : t - 1})

and

p (x_{N} | y_{1 : N})

obtained from Algorithm 1 and

p (y_{t : N} | x_{t})

obtained from Algorithm 2.

₂ Save the PDF

p (x_{N} | y_{1 : N})

₁₃ Output: The state estimation in (11), the covariance matrix of the estimation error

in (13), and the smoothing PDFs

p (x_{t} | y_{1 : N})

, for

t = 1, \dots, N

3.9. Computing the Smoothing Joint PDF $p (x_{t + 1}, x_{t} | y_{1 : N})$

In Theorems 2–4, the filtering and smoothing problems for quantized data are solved. However, many strategies for the identification of state-state space systems [8,9] require the joint PDF

p (x_{t + 1}, x_{t} | y_{1 : N})

, which for the quantized output data case is given by the following:

Theorem 5.

Consider the system in (4)–(6), the PDF

p (x_{t} | y_{1 : t})

and the backward prediction equation

p (y_{t + 1 : N} | x_{t + 1})

given in Theorems 2 and 3, respectively, and the PDF

p (x_{t + 1} | x_{t})

given in Equation (5). Then, for

t = N - 1, \dots, 1

, the joint PDF

p (x_{t + 1}, x_{t} | y_{1 : N})

is the GMM given by:

p (x_{t + 1}, x_{t} | y_{1 : N}) = \sum_{k = 1}^{S_{t, t + 1}} α_{k} N_{χ_{t}} ({\hat{χ}}_{t | N}^{k}, E_{t | N}^{k}),

(61)

where

S_{t, t + 1} = M_{t | t} S_{t + 1 | t + 1}

M_{t | t}

and

S_{t + 1 | t + 1}

are given in (33) and (42), respectively,

α_{k}

{\hat{χ}}_{t | N}^{k}

, and

E_{t | N}^{k}

are given in Appendix E, and

χ_{t}

is the extended vector given by:

χ_{t}^{T} = {[\begin{matrix} x_{t + 1}^{T} & x_{t}^{T} \end{matrix}]}^{T} .

(62)

At time

t = N

, the joint PDF

p (x_{N + 1}, x_{N} | y_{1 : N})

is given by (61) with:

\begin{matrix} S_{N + 1 | N + 1} = 1, ϵ_{N + 1 | N + 1}^{ℓ} = 1, λ_{N + 1 | N + 1}^{ℓ} = 1, \\ F_{N + 1 | N + 1}^{ℓ} = 0, G_{N + 1 | N + 1}^{ℓ T} = 0, H_{N + 1 | N + 1}^{ℓ} = 0 . \end{matrix}

(63)

Proof.

Using Bayes’ theorem, the joint PDF

p (x_{t + 1}, x_{t} | y_{1 : N})

can be obtained as:

p (x_{t + 1}, x_{t} | y_{1 : N}) = \frac{p (y_{t + 1 : N} | x_{t + 1}) p (x_{t + 1} | x_{t}) p (x_{t} | y_{1 : t})}{p (y_{t + 1 : N} | y_{1 : t})},

(64)

where

p (y_{t + 1 : N} | y_{1 : t})

is a normalization constant. The required PDFs

p (x_{t} | y_{1 : t})

and

p (x_{t + 1} | x_{t})

and backward-prediction equation

p (y_{t + 1 : N} | x_{t + 1})

are given in (33), (5) and (42), respectively. Using Lemma A2,

p (x_{t} | y_{1 : t})

in (33) can be rewritten as:

\begin{matrix} p (x_{t} | y_{1 : t}) = \sum_{τ = 1}^{M_{t | t}} \frac{γ_{t | t}^{τ}}{\sqrt{det \{2 π Σ_{t | t}^{τ}\}}} exp \{- \frac{1}{2} (x_{t}^{T} {(Σ_{t | t}^{τ})}^{- 1} x_{t} - 2 J_{t}^{τ T} x_{t} + L_{t}^{τ})\}, \end{matrix}

(65)

where

J_{t}^{τ T}

and

L_{t}^{τ}

are defined in (A44) and (A45), respectively. Considering

p (y_{t + 1 : N} | x_{t + 1})

p (x_{t} | y_{1 : t})

p (x_{t + 1} | x_{t})

, and the extended vector of the state given in (62), the argument of these three functions can be written as follows:

A_{1} = χ_{t}^{T} [\begin{matrix} 0 & 0 \\ 0 & {(Σ_{t | t}^{τ})}^{- 1} \end{matrix}] χ_{t} - 2 [\begin{matrix} 0 & J_{t}^{τ T} \end{matrix}] χ_{t} + L_{t}^{τ},

(66)

A_{2} = χ_{t}^{T} [\begin{matrix} F_{t + 1 | t + 1}^{ℓ} & 0 \\ 0 & 0 \end{matrix}] χ_{t} - 2 [\begin{matrix} G_{t + 1 | t + 1}^{ℓ T} & 0 \end{matrix}] χ_{t} + H_{t + 1 | t + 1}^{ℓ},

(67)

\begin{matrix} A_{3} & = χ_{t}^{T} [\begin{matrix} Q^{- 1} & - Q^{- 1} A \\ - A^{T} Q^{- 1} & A^{T} Q^{- 1} A \end{matrix}] χ_{t} \\ - 2 [\begin{matrix} u_{t}^{T} B^{T} Q^{- 1} & - u_{t}^{T} B^{T} Q^{- 1} A \end{matrix}] χ_{t} + u_{t}^{T} B^{T} Q^{- 1} B u_{t}, \end{matrix}

(68)

Then, adding these expressions, the following equation is derived:

\begin{matrix} p (x_{t + 1}, x_{t} | y_{1 : N}) & \propto \sum_{τ = 1}^{M_{t | t}} \sum_{ℓ = 1}^{S_{t + 1 | t + 1}} \frac{γ_{t | t}^{τ} ϵ_{t + 1 | t + 1}^{ℓ} λ_{t + 1 | t + 1}^{ℓ}}{{(2 π)}^{n} \sqrt{det \{Q\} det \{Σ_{t | t}^{τ}\}}} \\ \times exp \{- \frac{1}{2} (χ_{t}^{T} F_{t}^{ℓ τ} χ_{t} - 2 G_{t}^{ℓ τ T} χ_{t} + H_{t}^{ℓ τ})\}, \end{matrix}

(69)

where

F_{t}^{ℓ τ}

G_{t}^{ℓ τ T}

, and

H_{t}^{ℓ τ}

are defined in (A41), (A42), and (A43), respectively. Completing the square and expressing the double summation as a single one with a new index

k = (ℓ - 1) M_{t | t} + τ

, it is obtained:

\begin{matrix} p (x_{t + 1}, x_{t} | y_{1 : N}) \propto \sum_{k = 1}^{S_{t, t + 1}} {\bar{α}}_{k} N_{χ_{t}} ({\hat{χ}}_{t | N}^{k}, E_{t | N}^{k}), \end{matrix}

(70)

where

S_{t, t + 1} = M_{t | t} S_{t + 1 | t + 1}

{\bar{α}}_{k}

is defined in (A39), and

{\hat{χ}}_{t | N}^{k}

and

E_{t | N}^{k}

are defined in (A37) and (A38), respectively. Finally, the smoothing joint distribution in (64) for

t = N - 1, \dots, 1

is derived by computing the normalized wights as in (A36). Notice that, for

t = N

, the smoothing joint PDF

p (x_{N + 1}, x_{N} | y_{1 : N}) = p (x_{N + 1} | x_{N}) p (x_{N} | y_{1 : N})

is obtained as (61) considering the expressions given in (63). Then, the proof is finished. □

In Algorithm 4, the steps to implement the Gaussian sum smoothing for computing the joint PDF

p (x_{t + 1}, x_{t} | y_{1 : N})

are summarized.

Algorithm 4 Gaussian sum smoother to compute

p (x_{t + 1}, x_{t} | y_{1 : N})

for quantized output data

₁ Input: The PDF

p (x_{t} | y_{1 : t})

obtained in Algorithm 1,

p (y_{t + 1 : N} | x_{t + 1})

computed in

Algorithm 2, and the PDF

p (x_{t + 1} | x_{t})

given in Equation (5).

₁₁ Compute and store

p (x_{N + 1}, x_{N} | y_{1 : N})

according to Theorem 5 with

S_{N + 1 | N + 1}

ϵ_{N + 1 | N + 1}^{ℓ}

λ_{N + 1 | N + 1}^{ℓ}

F_{N + 1 | N + 1}^{ℓ}

G_{N + 1 | N + 1}^{ℓ T}

, and

H_{N + 1 | N + 1}^{ℓ}

given in (63).

₁₂ Output: The smoothing PDFs

p (x_{t + 1}, x_{t} | y_{1 : N})

, for

t = 1, \dots, N

4. Numerical Example

In this section, a numerical example to illustrate the benefits of this paper proposal is presented. Furthermore, a practical simulation example is studied: the problem of estimating the tank liquid level in a two-tank system is addressed. Typically, a numerical simulation approach is used for testing new algorithms and designs in applications of state estimation [74], control [75], and system identification [63], among others. This approach is used to evaluate the performance of the estimation algorithms, in order to avoid safety problems that can occur in a real-world processes.

To evaluate the performance of the proposed filtering and smoothing methods for quantized data, a comparison with three classical techniques is presented: standard Kalman filtering and smoothing [76], quantized Kalman filtering and smoothing [49,77], and particle filtering and smoothing [31]. Notice that it is also possible to implement a version of particle filtering and smoothing algorithms using the approximation of

p (y_{t} | x_{t})

given in (21). However, in the numerical examples run for filtering, the computation time was high. In order to validate the approximation in Theorem 1, the standard particle filtering and smoothing algorithms with a large number of particles were used, where

p (y_{t} | x_{t})

was computed using the integral in (20) from the MATLAB function mvncdf, which computes the multivariate normal cumulative distribution function. The true filtering and smoothing PDFs were considered to be provided by the particle filter and smoother with 20,000 particles, which was defined as the ground truth.

In the following examples, the discrete-time system in state-space form given in (1)–(3) is used with:

y_{t} = Δ_{q} round \{z_{t} / Δ_{q}\},

(71)

where the quantizer is defined in terms of the round function in MATLAB and the quantization step

Δ_{q}

. The infimum and supremum values of the sets

J_{i}

defined in (8),

q_{i - 1}

and

q_{i}

, can be calculated for the infinite-level quantizer as

q_{i - 1} = y_{t} - 0.5 Δ_{q}

and

q_{i} = y_{t} + 0.5 Δ_{q}

for

i = \dots, 1, 2 \dots, L, \dots

The experiments were carried out on a computer with the following specifications: processor: Intel(R) Core(TM) i5-8300H CPU @ 2.30GHz, RAM memory: 8.00 GB, operating system: Windows 10 with MATLAB 2020b.

4.1. Example 1: First-Order System

In this example, the following state-space system was considered:

\begin{matrix} (72) & x_{t + 1} & = 0.9 x_{t} + 1.0 u_{t} + w_{t}, \\ (73) & z_{t} & = 2.0 x_{t} + 0.5 u_{t} + v_{t}, \end{matrix}

where

w_{t} \sim N_{w_{t}} (0, 1)

and

v_{t} \sim N_{v_{t}} (0, 0.5)

. Furthermore, the input was considered to be drawn from

N_{u_{t}} (0, 2)

and

x_{1} \sim N_{x_{1}} (1, 1)

. In Figure 5, the filtering PDFs for some time instants are shown. The quantization step used in this example and the number of Gaussian components to approximate

p (y_{t} | x_{t})

were considered to be

Δ_{q} = 7

and

K = 10

, respectively. Furthermore,

S_{red} = K

Figure 5 shows that the filtering PDFs obtained with our proposal, the Gaussian sum filter, were that best fit to the ground truth. In contrast, the PDFs obtained with the Kalman filter, quantized Kalman filter, and particle filter were different from the ground truth. Notice that the results obtained using the particle filter with 100 particles were good at

t = 3, 9, 40

, whilst the resulting PDFs for

t = 25, 32, 59, 72, 82, 99

differed slightly more from the ground truth. However, the performance of particle filtering can be improved by increasing the number of particles, which produces an increment in the execution time.

To compare the accuracy of the state estimation, 100 Monte Carlo trials were run, and the mean-squared error between the true state and the estimation obtained by the Kalman filter, quantized Kalman filter, Gaussian sum filter, and particle filter is computed as follows:

MSE = \frac{1}{R} \sum_{k = 1}^{R} ‖ x_{t} - {\hat{x}}_{t | t} ‖_{2}^{2},

(74)

where

x_{t}

is the true state (which in a real-word system is not available, but in simulations, it can be used to analyze the performance of the estimation techniques),

{\hat{x}}_{t | t}

is the estimation of the state system, and

‖ \cdot ‖_{2}^{2}

denotes the squared Euclidean norm. In Figure 6, the box plots corresponding to 100 Monte Carlo runs for

Δ_{q} = \{3, 5, 7\}

are shown.

Figure 6 shows that the MSE between the true state and the estimation obtained with the Kalman filter and quantized Kalman filter increased fast as

Δ_{q}

increased. However, the MSE increased slowly for the state estimation obtained with the Gaussian sum filter and particle filter.

In Figure 7, the smoothing PDFs at time

t = 100

are shown. The quantization step was considered to be chosen from

Δ_{q} = \{3, 5, 7\}

, and the number of Gaussian components to approximate

p (y_{t} | x_{t})

was chosen from

K = \{6, 8, 10\}

. Furthermore, the number of Gaussian components kept after the reduction procedure was considered as

S_{red} = K

. In order to obtain the adequate number of particles to compare the particle smoother with the Gaussian sum smoother, 100 Monte Carlo simulations were carried out to obtain the number of particles that yielded smoothing PDFs that were as close to the true PDFs as the Gaussian sum smoother using

K = \{6, 8, 10\}

. For comparison purposes, the particle smoother execution time corresponding to the time required to implement the particle smoother algorithm using the number of particles that produces a similar result to the Gaussian sum smoother with K components is defined as Par(K). The

L_{2}

-norm of the difference between the true and the estimated PDFs as the measure of similarity was used:

‖ q - \hat{q} ‖_{2} = {[\sum_{k = 1}^{M} {| q - \hat{q} |}^{2}]}^{1 / 2},

(75)

where q represents the true PDF and

\hat{q}

represents the estimated PDF, which was chosen so that

‖ q - \hat{q} ‖_{2} < 1 \times 10^{- 6}

. The approximated number of particles (labeled in each PDF in Figure 7) was used to compare the execution time of both algorithms, the Gaussian sum smoother and particle smoother, and the results are shown in Table 2. Figure 7 shows that the smoothing PDFs obtained with our proposal using a small number of Gaussian components and the PDFs obtained using the particle smoother with a large number of particles fit the ground truth. In contrast, the PDFs obtained with the Kalman smoother and quantized Kalman smoother were different from the ground truth. Furthermore, the execution time in Table 2 shows that the required time to perform the Gaussian sum smoother was less compared to the time to perform the particle smoother, which needs a large number of particles to produce a result similar to the Gaussian sum smoother.

From the results shown in Figure 5, Figure 6 and Figure 7 and Table 2, it can be concluded that:

The filtering and smoothing PDFs are non-Gaussian, although the process and output noises in (1) and (2) are Gaussian distributed;
The accuracy of the standard and quantized Kalman filtering and smoothing decreased as the quantization step increased;
The state estimates obtained with particle filter and smoother were similar to the results obtained using the Gaussian sum filter and smoother. However, the characterization of the filtering and smoothing PDFs using the Gaussian sum filter and smoother were better than the PDF obtained by the particle filter and smoother. Notice that a correct characterization of a PDF is important when high-order moments need to be computed, especially in system identification tasks;
In order to implement the Gaussian sum filter and smoother, the parameters K (the number that defines the quality of the $p (y_{t} | x_{t})$ approximation) and $S_{red}$ (the Gaussian components kept after the Gaussian sum reduction algorithm) need to be chosen by the user. These parameters can be found in a simulation study by the trial and error approach and should be set by a trade-off between the time complexity and the accuracy of the estimation. A large value of K produces an accurate estimate, but a high computational load;
The larger the quantization step is, the larger the number of Gaussian components, K, needed to approximate $p (y_{t} | x_{t})$ in order to obtain an accurate approximation. However, for a large quantization step, the number K needed to obtain a good approximation of the filtering and smoothing PDFs is relatively small compared to the number of particles required to obtain similar results using the particle filter and smoother;
The maximum number of Gaussian components kept after the Gaussian reduction procedure is important for the accuracy of the approximation. In the simulations, $S_{red} = K$ was used. Furthermore, it was noticed that once an adequate $S_{red}$ was defined, incrementing this value did not produce a significant improvement in the estimation. However, this increment in $S_{red}$ was really critical for the resulting numerical complexity of the algorithm (and hence, the execution time), which increased since the Gaussian sum reduction procedure (e.g., Kullback–Leibler reduction) utilized more time to reduce a large amount of Gaussian components;
The Gaussian sum smoother execution time for all values of $Δ_{q}$ was small. This occurred because in each case, a relatively small number of Gaussian components to approximate $p (y_{t} | x_{t})$ were used. However, the particle smoother execution time is variable for different values of $Δ_{q}$ . As $Δ_{q}$ decreased, the $L_{2}$ -norm between the Gaussian sum smoother and the ground truth decreased, and a larger number of particles to obtain a comparable $L_{2}$ -norm between the particle smoother and the ground truth were required.

4.2. Real-World Application: Tank Liquid Level

In this example, the problem of estimating the tank liquid level in a two-tank system by using the measurements taken by a low-cost sensor based on a variable resistor that is attached to an arm with a floater is considered; see Figure 8.

A linearized model of this system can be found in [78]. Here, it was assumed that

h_{2}

can be measured and

h_{1}

cannot. The discrete-time version of the model in [78] with sample time

0.1

s was considered:

\begin{matrix} (76) & x_{t + 1} & = [\begin{matrix} 0.9959 & 0.0041 \\ 0.0041 & 0.9959 \end{matrix}] x_{t} + [\begin{matrix} 0.0998 & - 0.0002 \\ 0.0002 & - 0.0998 \end{matrix}] u_{t} + w_{t}, \\ (77) & z_{t} & = [\begin{matrix} 0 & 1.0 \end{matrix}] x_{t} + v_{t}, \end{matrix}

where

x_{t} = {[h_{1} h_{2}]}^{T}

and

u_{t} = {[f_{1} - k f_{2} + k]}^{T}

with

k = 0.4111

. The sensor measures

h_{2}

vary its resistance in discrete steps, with minimum and maximum values

β_{1} = 2

and

β_{L} = 10

. To simulate this system, it was considered that:

w_{t} \sim N_{w_{t}} (0, 0.001 I_{2})

v_{t} \sim N_{v_{t}} (0, 0.0001)

; the input

{[f_{1} f_{2}]}^{T}

was drawn from

N ({[10 2]}^{T}, 10 I_{2})

; the initial condition

x_{1} \sim N_{x_{1}} ({[10 5]}^{T}, 0.01 I_{2})

. For this example, 100 Monte Carlo runs were simulated. In Figure 9 (left), the output

z_{t}

that corresponds to the values of

h_{2}

that are nonquantized and the output

y_{t}

that corresponds to the measurements given by the sensor (for one of the Monte Carlo runs) are shown. In this figure, the level of quantization in the measurements can be observed. In Figure 9 (right), the MSE between the true and estimated state is shown. It was observed that the proposal presented in this paper, the Gaussian sum filter, yielded the most accurate estimation of

h_{1}

, followed by the particle filter, Kalman filter, and quantized Kalman filter. In this example,

K = 50

and 1000 particles were used to implement the Gaussian sum filter and particle filter, respectively.

In this example, a relatively high number of Gaussian components

(K = 50)

were required to obtain a good estimation of the filtering and smoothing distributions, and hence of the state. This produced an increment in the execution time since the Gaussian sum reduction algorithm needed more time to deal with a high number of Gaussian components in every iteration. This resulted in similar execution times for our proposed algorithm and the traditional particle filter. However, the execution time of the Gaussian sum filter and smoother was smaller than the execution time of the ground truth (particle filter and smoother with 20,000 particles).

5. Conclusions

In this paper, Gaussian sum filtering and smoothing algorithms for linear-time-invariant state-space systems with quantized output data were developed. An approximation of the integral equation that defines the probability mass function of the quantized data given the current state,

p (y_{t} | x_{t})

, using Gaussian quadrature was considered. This approximation naturally yielded an explicit mathematical model with a GMM structure for this probability function. Using the approximation of

p (y_{t} | x_{t})

summarized in Theorem 1, it was possible to solve in closed form the general equations of filtering and smoothing to deal with quantized data. This fact allowed for a closed-form expression of the system state estimators given the quantized data

{\hat{x}}_{t | t}

and

{\hat{x}}_{t | N}

. Via numerical simulations, it was shown that approximating

p (y_{t} | x_{t})

with a small number of Gaussian components was adequate, yielding an approximation comparable to the true filtering and smoothing PDFs given by the particle approach (using a large number of particles, namely 20,000 particles). This reduced number of Gaussian components allowed for a low computational load, especially when the system order increased. In addition, our results showed overall less computational load for our proposed techniques since the number of Gaussian components was considerably less than the number of particles used in particle filtering and smoothing.

The proposed Gaussian sum filtering and smoothing algorithms can be utilized, in principle, to develop system identification algorithms and control strategies having quantized output measurements.

Author Contributions

Conceptualization, A.L.C., R.A. and J.C.A.; methodology, A.L.C. and J.C.A.; software, A.L.C.; validation, R.C. and B.I.G.; formal analysis, A.L.C., R.A., R.C., B.I.G. and J.C.A.; investigation, A.L.C., R.A. and J.C.A.; resources, J.C.A.; writing—original draft preparation, A.L.C., R.A., R.C., B.I.G. and J.C.A.; writing—review and editing, A.L.C., R.C., B.I.G. and J.C.A.; visualization, R.C. and B.I.G.; supervision, J.C.A. All authors have read and agreed to the published version of the manuscript.

Funding

PIIC program of DGP at Universidad Técnica Federico Santa María No. 062/2018 and 035/2021. Grants ANID-Fondecyt 1211630 and 11201187, ANID-Basal Project FB0008 (AC3E). Chilean National Agency for Research and Development (ANID) Scholarship Program/Doctorado Nacional/2020-21202410.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

PDF/PMF:	Probability density/mass function
GMM:	Gaussian mixture model
GSF/GSS:	Gaussian sum filter/smoother
PF/PS:	Particle filter/smoother
KF/KS:	Kalman filter/smoother
QKF/QKS:	Quantized Kalman filter/smoother
FLQ/ILQ:	Finite-/infinite-level quantizer
MSE:	Mean squared error
$S_{red}$ :	Number of Gaussian components kept after the reduction procedure
GT:	Ground truth obtain by using the particle filter/smoother with a large number
	of particles
K:	Order of the polynomials in the Gauss–Legendre quadrature

Appendix A. Technical Lemmata

Lemma A1.

The PDF

N_{y} (C x + μ, R)

of the random variable

Y \in R^{p}

can be rewritten as follows:

N_{y} (C x + μ, R) = \frac{1}{\sqrt{det \{2 π R\}}} exp \{- \frac{1}{2} (x^{T} F x - 2 G^{T} x + H)\},

(A1)

where

F = C^{T} R^{- 1} C

G^{T} = {(y - μ)}^{T} R^{- 1} C

, and

H = {(y - μ)}^{T} R^{- 1} (y - μ)

Proof.

Directly expand the exponential argument and reorder the terms in the variable x. □

Lemma A2.

The PDF

N_{x} (A w + ν, Q)

of the random variable

X \in R^{n}

can be rewritten as follows:

N_{x} (A w + ν, Q) = \frac{1}{\sqrt{det \{2 π Q\}}} exp \{- \frac{1}{2} (x^{T} Q^{- 1} x - 2 J^{T} x + L)\},

(A2)

where

J^{T} = {(A w + ν)}^{T} Q^{- 1}

and

L = {(A w + ν)}^{T} Q^{- 1} (A w + ν)

Proof.

Directly expand the exponential argument and reorder the terms in the variable x. □

Lemma A3.

Consider the backward filter function

p (y_{t : N} | x_{t})

given in (42). Its GMM structure is given by:

p (y_{t : N} | x_{t}) = \sum_{k = 1}^{S_{t | t}} σ_{t | t}^{k} N_{x_{t}} (ν_{t | t}^{k}, Ω_{t | t}^{k}),

(A3)

where

σ_{t | t}^{k} = {\bar{σ}}_{t | t}^{k} / \sum_{s = 1}^{S_{t | t}} {\bar{σ}}_{t | t}^{s}

is the normalized mixing weight and

ν_{t | t}^{k} = {(F_{t | t}^{k})}^{- 1} G_{t | t}^{k}

and

Ω_{t | t}^{k} = {(F_{t | t}^{k})}^{- 1},

are the mean vector and the covariance matrix, respectively, with:

{\bar{σ}}_{t | t}^{k} = \frac{ϵ_{t | t}^{k} λ_{t | t}^{k} exp \{- \frac{1}{2} (H_{t | t}^{k} - G_{t | t}^{k T} {(F_{t | t}^{k})}^{- 1} G_{t | t}^{k})\}}{{(2 π)}^{- \frac{n}{2}} \sqrt{det \{F_{t | t}^{k}\}}} .

(A4)

On the other hand, consider the GMM structure given by (54). Then, the backward filter form in (42) is obtained using the following:

\begin{matrix} S_{t | t} & = S_{red}, \\ F_{t | t}^{k} & = {(U_{t | t}^{k})}^{- 1}, \end{matrix} \begin{matrix} ϵ_{t | t}^{k} & = δ_{t | t}^{k}, \\ G_{t | t}^{k T} & = {(z_{t | t}^{k})}^{T} {(U_{t | t}^{k})}^{- 1}, \end{matrix} \begin{matrix} λ_{t | t}^{k} & = {(det \{2 π U_{t | t}^{k}\})}^{- 1 / 2}, \\ H_{t | t}^{k} & = {(z_{t | t}^{k})}^{T} {(U_{t | t}^{k})}^{- 1} (z_{t | t}^{k}) . \end{matrix}

(A5)

Proof.

Directly by using Lemma A1. □

Appendix B. Quantities to Perform the Gaussian Sum Filter

For each two-tuple

(τ, ℓ)

, where

τ = 1, \dots, K

and

ℓ = 1, \dots, M_{t | t - 1}

, let k be a new index, so that

k = (ℓ - 1) K + τ

. Then:

\begin{matrix} (A6) & γ_{t | t}^{k} & = \frac{{\bar{γ}}_{t | t}^{k}}{\sum_{s = 1}^{M_{t | t}} {\bar{γ}}_{t | t}^{s}}, \\ (A7) & {\hat{x}}_{t | t}^{k} & = {\hat{x}}_{t | t - 1}^{ℓ} + K_{t}^{ℓ} (η_{t}^{τ} - κ_{t}^{ℓ τ}), \\ (A8) & Σ_{t | t}^{k} & = (I - K_{t}^{ℓ} C) Σ_{t | t - 1}^{ℓ}, \end{matrix}

are the weights, means, and covariance matrices of the measurement-update equation in the Gaussian sum filter, with the following definitions:

\begin{matrix} (A9) & {\bar{γ}}_{t | t}^{k} & = ς_{t}^{τ} γ_{t | t - 1}^{ℓ} N_{η_{t}^{τ}} (κ_{t}^{ℓ τ}, V_{t}^{ℓ}), \\ (A10) & K_{t}^{ℓ} & = Σ_{t | t - 1}^{ℓ} C^{T} {(V_{t}^{ℓ})}^{- 1}, \\ (A11) & κ_{t}^{ℓ τ} & = C {\hat{x}}_{t | t - 1}^{ℓ} + D u_{t} + μ_{t}^{τ}, \\ (A12) & V_{t}^{ℓ} & = R + C Σ_{t | t - 1}^{ℓ} C^{T}, \end{matrix}

where

ς_{t}^{τ}

η_{t}^{τ}

, and

μ_{t}^{τ}

are defined in Table 1, and:

\begin{matrix} (A13) & γ_{t + 1 | t}^{k} & = γ_{t | t}^{k}, \\ (A14) & {\hat{x}}_{t + 1 | t}^{k} & = A {\hat{x}}_{t | t}^{k} + B u_{t}, \\ (A15) & Σ_{t + 1 | t}^{k} & = Q + A Σ_{t | t}^{k} A^{T}, \end{matrix}

are the weights, means, and covariance matrices of the time-update equation in the Gaussian sum filter.

Appendix C. Quantities to Perform the Backward Filter

For each two-tuple

(τ, ℓ)

, where

τ = 1, \dots, K

and

ℓ = 1, \dots, S_{t | t + 1}

, let k be a new index, so that

k = (ℓ - 1) K + τ

. Then:

\begin{matrix} (A16) & ϵ_{t | t}^{k} & = ς_{t}^{τ} ϵ_{t | t + 1}^{ℓ}, \\ (A17) & λ_{t | t}^{k} & = {(det \{2 π R\})}^{- 1 / 2} λ_{t | t + 1}^{ℓ}, \\ (A18) & θ_{t}^{τ} & = η_{t}^{τ} - D u_{t} - μ_{t}^{τ}, \\ (A19) & F_{t | t}^{k} & = F_{t | t + 1}^{ℓ} + C^{T} R^{- 1} C, \\ (A20) & G_{t | t}^{k T} & = G_{t | t + 1}^{ℓ T} + θ_{t}^{τ T} R^{- 1} C, \\ (A21) & H_{t | t}^{k} & = H_{t | t + 1}^{ℓ} + θ_{t}^{τ T} R^{- 1} θ_{t}^{τ} . \end{matrix}

are the quantities needed to compute the backward-measurement-update equation in the backward filter, where

ς_{t}^{τ}

η_{t}^{τ}

, and

μ_{t}^{τ}

are defined in Table 1, and:

\begin{matrix} (A22) & ϵ_{t | t + 1}^{k} & = ϵ_{t + 1 | t + 1}^{k}, \\ (A23) & λ_{t | t + 1}^{k} & = {(det \{Q\} det \{F_{q k}\})}^{- 1 / 2} λ_{t + 1 | t + 1}^{k}, \\ (A24) & F_{t | t + 1}^{k} & = A^{T} M_{q k} A, \\ (A25) & G_{t | t + 1}^{k T} & = G_{t + 1 | t + 1}^{k T} F_{q k}^{- 1} Q^{- 1} A - u_{t}^{T} B^{T} M_{q k} A, \\ (A26) & H_{t | t + 1}^{k} & = H_{t + 1 | t + 1}^{k} - G_{t + 1 | t + 1}^{k T} F_{q k}^{- 1} G_{t + 1 | t + 1}^{k} + u_{t}^{T} B^{T} M_{q k} B u_{t} - 2 u_{t}^{T} B^{T} Q^{- 1} F_{q k}^{- 1} G_{t + 1 | t + 1}^{k}, \end{matrix}

are the quantities needed to compute the backward prediction equation in the backward filter, where

F_{q k} = F_{t + 1 | t + 1}^{k} + Q^{- 1}

and

M_{q k} = Q^{- 1} - Q^{- 1} F_{q k}^{- 1} Q^{- 1}

Appendix D. Quantities to Perform the Gaussian Sum Smoother

For each two-tuple

(τ, ℓ)

, where

τ = 1, \dots, M_{t | t - 1}

and

ℓ = 1, \dots, S_{red}

, let k be a new index, so that

k = (ℓ - 1) M_{t | t - 1} + τ

. Then:

\begin{matrix} (A27) & ϵ_{t | N}^{k} & = \frac{{\bar{ϵ}}_{t | N}^{k}}{\sum_{s = 1}^{S_{t | N}} {\bar{ϵ}}_{t | N}^{s}}, \\ (A28) & {\hat{x}}_{t | N}^{k} & = {(L_{t}^{ℓ τ})}^{- 1} ρ_{t}^{ℓ τ}, \\ (A29) & Σ_{t | N}^{k} & = {(L_{t}^{ℓ τ})}^{- 1}, \end{matrix}

are the weights, means, and covariance matrices of the smoothing PDF

p (x_{t} | y_{1 : N})

, where:

\begin{matrix} (A30) & {\bar{ϵ}}_{t | N}^{k} & = \frac{γ_{t | t - 1}^{τ} δ_{t | t}^{ℓ} exp \{- \frac{1}{2} [ϕ_{1 t}^{ℓ} + ϕ_{2 t}^{τ} - ϕ_{3 t}^{ℓ τ}]\}}{{(2 π)}^{\frac{n}{2}} \sqrt{det \{L_{t}^{ℓ τ}\} det \{U_{t | t}^{ℓ}\} det \{Σ_{t | t - 1}^{τ}\}}}, \\ (A31) & ρ_{t}^{ℓ τ} & = [{(U_{t | t}^{ℓ})}^{- 1} z_{t | t}^{ℓ} + {(Σ_{t | t - 1}^{τ})}^{- 1} {\hat{x}}_{t | t - 1}^{τ}], \\ (A32) & L_{t}^{ℓ τ} & = [{(U_{t | t}^{ℓ})}^{- 1} + {(Σ_{t | t - 1}^{τ})}^{- 1}], \\ (A33) & ϕ_{3 t}^{ℓ τ} & = {(ρ_{t}^{ℓ τ})}^{T} {(L_{t}^{ℓ τ})}^{- 1} (ρ_{t}^{ℓ τ}), \\ (A34) & ϕ_{1 t}^{ℓ} & = {(z_{t | t}^{ℓ})}^{T} {(U_{t | t}^{ℓ})}^{- 1} (z_{t | t}^{ℓ}), \\ (A35) & ϕ_{2 t}^{τ} & = {({\hat{x}}_{t | t - 1}^{τ})}^{T} {(Σ_{t | t - 1}^{τ})}^{- 1} ({\hat{x}}_{t | t - 1}^{τ}), \end{matrix}

γ_{t | t - 1}^{τ}

{\hat{x}}_{t | t - 1}^{τ}

and

Σ_{t | t - 1}^{τ}

are obtained from the time-update step of Theorem 2 (34), and

δ_{t | t}^{ℓ}

z_{t | t}^{ℓ}

, and

U_{t | t}^{ℓ}

are obtained from the reduced measurement-update step of the backward-filtering algorithm in (54).

Appendix E. Quantities to Compute p(x_t+1, x_t|y_1:N)

For each two-tuple

(τ, ℓ)

, where

τ = 1, \dots, M_{t | t}

and

ℓ = 1, \dots, S_{t + 1 | t + 1}

, let k be a new index, so that

k = (ℓ - 1) M_{t | t} + τ

. Then:

\begin{matrix} (A36) & α_{k} & = \frac{{\bar{α}}_{k}}{\sum_{s = 1}^{S_{t, t + 1}} {\bar{α}}_{s}}, \\ (A37) & {\hat{χ}}_{t | N}^{k} & = {(F_{t}^{ℓ τ})}^{- 1} G_{t}^{ℓ τ}, \\ (A38) & E_{t | N}^{k} & = {(F_{t}^{ℓ τ})}^{- 1}, \end{matrix}

are the weights, means, and covariance matrices of the joint PDF

p (x_{t + 1}, x_{t} | y_{1 : N})

, where:

\begin{matrix} (A39) & {\bar{α}}_{k} & = \frac{γ_{t | t}^{τ} ϵ_{t + 1 | t + 1}^{ℓ} λ_{t + 1 | t + 1}^{ℓ} S_{t}^{ℓ τ}}{\sqrt{det \{Q\} det \{F_{t}^{ℓ τ}\} det \{Σ_{t | t}^{τ}\}}}, \\ (A40) & S_{t}^{ℓ τ} & = exp \{- \frac{1}{2} (H_{t}^{ℓ τ} - G_{t}^{ℓ τ T} {(F_{t}^{ℓ τ})}^{- 1} G_{t}^{ℓ τ})\}, \\ (A41) & F_{t}^{ℓ τ} & = [\begin{matrix} Q^{- 1} + F_{t + 1 | t + 1}^{ℓ} & - Q^{- 1} A \\ - A^{T} Q^{- 1} & {(Σ_{t | t}^{τ})}^{- 1} + A^{T} Q^{- 1} A \end{matrix}], \\ (A42) & G_{t}^{ℓ τ T} & = [\begin{matrix} G_{t + 1 | t + 1}^{ℓ T} + u_{t}^{T} B^{T} Q^{- 1} & J_{t}^{τ T} - u_{t}^{T} B^{T} Q^{- 1} A \end{matrix}], \\ (A43) & H_{t}^{ℓ τ} & = H_{t + 1 | t + 1}^{ℓ} + L_{t}^{τ} + u_{t}^{T} B^{T} Q^{- 1} B u_{t}, \end{matrix}

with:

\begin{matrix} (A44) & J_{t}^{τ T} & = {({\hat{x}}_{t | t}^{τ})}^{T} {(Σ_{t | t}^{τ})}^{- 1}, \\ (A45) & L_{t}^{τ} & = {({\hat{x}}_{t | t}^{τ})}^{T} {(Σ_{t | t}^{τ})}^{- 1} ({\hat{x}}_{t | t}^{τ}), \end{matrix}

γ_{t | t}^{τ}

{\hat{x}}_{t | t}^{τ}

, and

Σ_{t | t}^{τ}

are obtained from the measurement-update step in Theorem 2, and

ϵ_{t + 1 | t + 1}^{ℓ}

λ_{t + 1 | t + 1}^{ℓ}

F_{t + 1 | t + 1}^{ℓ}

G_{t + 1 | t + 1}^{ℓ T}

, and

H_{t + 1 | t + 1}^{ℓ}

are obtained from the reduced backward-measurement-update step in Theorem 3.

References

Kalman, R.E. A New Approach to Linear Filtering and Prediction Problems. Trans. ASME-J. Basic Eng. 1960, 82, 35–45. [Google Scholar] [CrossRef] [Green Version]
Kamen, E.W.; Su, J.K. Introduction to Optimal Estimation; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1999. [Google Scholar]
Anderson, B.D.O.; Moore, J.B. Optimal Filtering; Prentice-Hall, Inc.: Hoboken, NJ, USA, 1979. [Google Scholar]
Leong, A.S.; Dey, S.; Quevedo, D.E. Transmission scheduling for remote state estimation and control with an energy harvesting sensor. Automatica 2018, 91, 54–60. [Google Scholar] [CrossRef]
Liang, H.; Guo, X.; Pan, Y.; Huang, T. Event-Triggered Fuzzy Bipartite Tracking Control for Network Systems Based on Distributed Reduced-Order Observers. IEEE Trans. Fuzzy Syst. 2021, 29, 1601–1614. [Google Scholar] [CrossRef]
Liu, L.; Gao, T.; Liu, Y.J.; Tong, S.; Chen, C.L.P.; Ma, L. Time-varying IBLFs-based adaptive control of uncertain nonlinear systems with full state constraints. Automatica 2021, 129, 109595. [Google Scholar] [CrossRef]
Liu, L.; Liu, Y.J.; Chen, A.; Tong, S.; Chen, C.L.P. Integral Barrier Lyapunov function-based adaptive control for switched nonlinear systems. Sci. China Inf. Sci. 2020, 63, 132203. [Google Scholar] [CrossRef] [Green Version]
Gibson, S.; Ninness, B. Robust maximum-likelihood estimation of multivariable dynamic systems. Automatica 2005, 41, 1667–1682. [Google Scholar] [CrossRef]
Agüero, J.C.; Tang, W.; Yuz, J.I.; Delgado, R.; Goodwin, G.C. Dual time–frequency domain system identification. Automatica 2012, 48, 3031–3041. [Google Scholar] [CrossRef]
Kaltiokallio, O.; Hostettler, R.; Yiğitler, H.; Valkama, M. Unsupervised Learning in RSS-Based DFLT Using an EM Algorithm. Sensors 2021, 21, 5549. [Google Scholar] [CrossRef] [PubMed]
Zhao, J.; Netto, M.; Huang, Z.; Yu, S.S.; Gómez-Expósito, A.; Wang, S.; Kamwa, I.; Akhlaghi, S.; Mili, L.; Terzija, V.; et al. Roles of Dynamic State Estimation in Power System Modeling, Monitoring and Operation. IEEE Trans. Power Syst. 2021, 36, 2462–2472. [Google Scholar] [CrossRef]
Ji, X.; Yin, Z.; Zhang, Y.; Wang, M.; Zhang, X.; Zhang, C.; Wang, D. Real-time robust forecasting-aided state estimation of power system based on data-driven models. Int. J. Electr. Power Energy Syst. 2021, 125, 106412. [Google Scholar] [CrossRef]
Bonvini, M.; Sohn, M.D.; Granderson, J.; Wetter, M.; Piette, M.A. Robust on-line fault detection diagnosis for HVAC components based on nonlinear state estimation techniques. Appl. Energy 2014, 124, 156–166. [Google Scholar] [CrossRef]
Nemati, F.; Safavi Hamami, S.M.; Zemouche, A. A nonlinear observer-based approach to fault detection, isolation and estimation for satellite formation flight application. Automatica 2019, 107, 474–482. [Google Scholar] [CrossRef]
Jeong, H.; Park, B.; Park, S.; Min, H.; Lee, S. Fault detection and identification method using observer-based residuals. Reliab. Eng. Syst. Saf. 2019, 184, 27–40. [Google Scholar] [CrossRef]
Noshad, Z.; Javaid, N.; Saba, T.; Wadud, Z.; Saleem, M.Q.; Alzahrani, M.E.; Sheta, O.E. Fault Detection in Wireless Sensor Networks through the Random Forest Classifier. Sensors 2019, 19, 1568. [Google Scholar] [CrossRef] [Green Version]
Huang, C.; Shen, B.; Zou, L.; Shen, Y. Event-Triggering State and Fault Estimation for a Class of Nonlinear Systems Subject to Sensor Saturations. Sensors 2021, 21, 1242. [Google Scholar] [CrossRef]
Corbetta, M.; Sbarufatti, C.; Giglio, M.; Todd, M.D. Optimization of nonlinear, non-Gaussian Bayesian filtering for diagnosis and prognosis of monotonic degradation processes. Mech. Syst. Signal Process. 2018, 104, 305–322. [Google Scholar] [CrossRef]
Ahwiadi, M.; Wang, W. An Adaptive Particle Filter Technique for System State Estimation and Prognosis. IEEE Trans. Instrum. Meas. 2020, 69, 6756–6765. [Google Scholar] [CrossRef]
Ding, D.; Han, Q.L.; Ge, X.; Wang, J. Secure State Estimation and Control of Cyber-Physical Systems: A Survey. IEEE Trans. Syst. Man Cybern. Syst. 2021, 51, 176–190. [Google Scholar] [CrossRef]
McLaughlin, D. An integrated approach to hydrologic data assimilation: Interpolation, smoothing, and filtering. Adv. Water Resour. 2002, 25, 1275–1286. [Google Scholar] [CrossRef]
Cosme, E.; Verron, J.; Brasseur, P.; Blum, J.; Auroux, D. Smoothing problems in a Bayesian framework and their linear Gaussian solutions. Mon. Weather Rev. 2012, 140, 683–695. [Google Scholar] [CrossRef]
Stone, L.D.; Streit, R.L.; Corwin, T.L.; Bell, K.L. Bayesian Multiple Target Tracking; Artech House: Norwood, MA, USA, 2013. [Google Scholar]
Schizas, I.D.; Giannakis, G.B.; Roumeliotis, S.I.; Ribeiro, A. Consensus in Ad Hoc WSNs With Noisy Links—Part II: Distributed Estimation and Smoothing of Random Signals. IEEE Trans. Signal Process. 2008, 56, 1650–1666. [Google Scholar] [CrossRef]
Liu, J.; Liu, Y.; Dong, K.; Ding, Z.; He, Y. A Novel Distributed State Estimation Algorithm with Consensus Strategy. Sensors 2019, 19, 2134. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Han, Y.; Cui, M.; Liu, S. Optimal Sensor and Relay Nodes Power Scheduling for Remote State Estimation with Energy Constraint. Sensors 2020, 20, 1073. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Chiang, K.W.; Duong, T.T.; Liao, J.K.; Lai, Y.C.; Chang, C.C.; Cai, J.M.; Huang, S.C. On-Line Smoothing for an Integrated Navigation System with Low-Cost MEMS Inertial Sensors. Sensors 2012, 12, 17372–17389. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Rostami Shahrbabaki, M.; Safavi, A.A.; Papageorgiou, M.; Papamichail, I. A data fusion approach for real-time traffic state estimation in urban signalized links. Transp. Res. Part C Emerg. Technol. 2018, 92, 525–548. [Google Scholar] [CrossRef]
Ahmed, A.; Naqvi, S.A.A.; Watling, D.; Ngoduy, D. Real-Time Dynamic Traffic Control Based on Traffic-State Estimation. Transp. Res. Rec. 2019, 2673, 584–595. [Google Scholar] [CrossRef]
Schreiter, T.; van Lint, H.; Treiber, M.; Hoogendoorn, S. Two fast implementations of the Adaptive Smoothing Method used in highway traffic state estimation. In Proceedings of the 13th International IEEE Conference on Intelligent Transportation Systems, Funchal, Portugal, 19–22 September2010; pp. 1202–1208. [Google Scholar]
Särkkä, S. Bayesian Filtering and Smoothing; Cambridge University Press: Cambridge, UK, 2013; Volume 3. [Google Scholar]
Widrow, B.; Kollár, I. Quantization Noise: Roundoff Error in Digital Computation, Signal Processing, Control, and Communications; Cambridge University Press: Cambridge, UK, 2008. [Google Scholar]
Msechu, E.J.; Giannakis, G.B. Sensor-Centric Data Reduction for Estimation With WSNs via Censoring and Quantization. IEEE Trans. Signal Process. 2012, 60, 400–414. [Google Scholar] [CrossRef]
Wang, L.Y.; Yin, G.G.; Zhang, J. System identification using binary sensors. IEEE Trans. Automat. Contr. 2003, 48, 1892–1907. [Google Scholar] [CrossRef] [Green Version]
Curry, R.E. Estimation and Control with Quantized Measurements; MIT Press: Cambridge, MA, USA, 1970. [Google Scholar]
Liu, S.; Wang, Z.; Hu, J.; Wei, G. Protocol-based extended Kalman filtering with quantization effects: The Round-Robin case. Int. J. Robust Nonlinear Control 2020, 30, 7927–7946. [Google Scholar] [CrossRef]
Malyavej, V.; Savkin, A.V. The problem of optimal robust Kalman state estimation via limited capacity digital communication channels. Syst. Control Lett. 2005, 54, 283–292. [Google Scholar] [CrossRef]
Farhadi, A.; Charalambous, C.D. Stability and reliable data reconstruction of uncertain dynamic systems over finite capacity channels. Automatica 2010, 46, 889–896. [Google Scholar] [CrossRef]
Silva, E.I.; Agüero, J.C.; Goodwin, G.C.; Lau, K.; Wang, M. The SNR Approach to Networked Control. In The Control Handbook: Control System Applications; Levine, W.S., Ed.; CRC Press: Boca Raton, FL, USA, 2011; Chapter 25; pp. 1–26. [Google Scholar]
Godoy, B.I.; Agüero, J.C.; Carvajal, R.; Goodwin, G.C.; Yuz, J.I. Identification of sparse FIR systems using a general quantisation scheme. Int. J. Control 2014, 87, 874–886. [Google Scholar] [CrossRef]
Li, Z.M.; Chang, X.H.; Yu, L. Robust quantized H∞ filtering for discrete-time uncertain systems with packet dropouts. Appl. Math. Comput. 2016, 275, 361–371. [Google Scholar] [CrossRef]
Li, S.; Sauter, D.; Xu, B. Fault isolation filter for networked control system with event-triggered sampling scheme. Sensors 2011, 11, 557–572. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Zhang, L.; Liang, H.; Sun, Y.; Ahn, C.K. Adaptive Event-Triggered Fault Detection Scheme for Semi-Markovian Jump Systems With Output Quantization. IEEE Trans. Syst. Man, Cybern. Syst. 2021, 51, 2370–2381. [Google Scholar] [CrossRef]
Zhang, X.; Han, Q.; Ge, X.; Ding, D.; Ding, L.; Yue, D.; Peng, C. Networked control systems: A survey of trends and techniques. IEEE/CAA J. Autom. Sin. 2020, 7, 1–17. [Google Scholar] [CrossRef]
Goodwin, G.C.; Haimovich, H.; Quevedo, D.E.; Welsh, J.S. A moving horizon approach to Networked Control system design. IEEE Trans. Autom. Control 2004, 49, 1427–1445. [Google Scholar] [CrossRef]
Gustafsson, F.; Karlsson, R. Statistical results for system identification based on quantized observations. Automatica 2009, 45, 2794–2801. [Google Scholar] [CrossRef] [Green Version]
Wang, L.Y.; Yin, G.G.; Zhang, J.; Zhao, Y. System Identification with Quantized Observations; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
Marelli, D.E.; Godoy, B.I.; Goodwin, G.C. A scenario-based approach to parameter estimation in state-space models having quantized output data. In Proceedings of the 49th IEEE Conference on Decision and Control (CDC), Atlanta, GA, USA, 15–17 December 2010; pp. 2011–2016. [Google Scholar]
Leong, A.S.; Dey, S.; Nair, G.N. Quantized Filtering Schemes for Multi-Sensor Linear State Estimation: Stability and Performance Under High Rate Quantization. IEEE Trans. Signal Process. 2013, 61, 3852–3865. [Google Scholar] [CrossRef]
Li, D.; Kar, S.; Alsaadi, F.E.; Dobaie, A.M.; Cui, S. Distributed Kalman Filtering with Quantized Sensing State. IEEE Trans. Signal Process. 2015, 63, 5180–5193. [Google Scholar] [CrossRef]
Rana, M.M.; Li, L. An Overview of Distributed Microgrid State Estimation and Control for Smart Grids. Sensors 2015, 15, 4302–4325. [Google Scholar] [CrossRef] [Green Version]
Chang, X.H.; Liu, Y. Robust H∞ Filtering for Vehicle Sideslip Angle With Quantization and Data Dropouts. IEEE Trans. Veh. Technol. 2020, 69, 10435–10445. [Google Scholar] [CrossRef]
Curiac, D. Towards wireless sensor, actuator and robot networks: Conceptual framework, challenges and perspectives. J. Netw. Comput. Appl. 2016, 63, 14–23. [Google Scholar] [CrossRef]
Allik, B.; Piovoso, M.J.; Zurakowski, R. Recursive estimation with quantized and censored measurements. In Proceedings of the 2016 American Control Conference (ACC), Boston, MA, USA, 6–8 July 2016; pp. 5130–5135. [Google Scholar]
Zhou, Y.; Li, J.; Wang, D. Unscented Kalman Filtering based quantized innovation fusion for target tracking in WSN with feedback. In Proceedings of the 2009 International Conference on Machine Learning and Cybernetics, Baoding, China, 12–15 July 2009; Volume 3, pp. 1457–1463. [Google Scholar]
Wigren, T. Approximate gradients, convergence and positive realness in recursive identification of a class of non-linear systems. Int. J. Adapt. Control Signal Process. 1995, 9, 325–354. [Google Scholar] [CrossRef]
Gustafsson, F.; Karlsson, R. Estimation based on Quantized Observations. IFAC Proc. Vol. 2009, 42, 78–83. [Google Scholar] [CrossRef]
Cedeño, A.L.; Albornoz, R.; Carvajal, R.; Godoy, B.I.; Agüero, J.C. On Filtering Methods for State-Space Systems having Binary Output Measurements. IFAC-PapersOnLine 2021, 54, 815–820. [Google Scholar] [CrossRef]
Gersho, A.; Gray, R.M. Vector Quantization and Signal Compression; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012; Volume 159. [Google Scholar]
Frühwirth, S.; Celeux, G.; Robert, C.P. Handbook of Mixture Analysis; CRC Press: Boca Raton, FL, USA, 2019. [Google Scholar]
Lo, J. Finite-dimensional sensor orbits and optimal nonlinear filtering. IEEE Trans. Inf. Theory 1972, 18, 583–588. [Google Scholar] [CrossRef]
Cedeño, A.L.; Orellana, R.; Carvajal, R.; Agüero, J.C. EM-based identification of static errors-in-variables systems utilizing Gaussian Mixture models. IFAC-PapersOnLine 2020, 53, 863–868. [Google Scholar] [CrossRef]
Orellana, R.; Carvajal, R.; Escárate, P.; Agüero, J.C. On the Uncertainty Identification for Linear Dynamic Systems Using Stochastic Embedding Approach with Gaussian Mixture Models. Sensors 2021, 21, 3837. [Google Scholar] [CrossRef] [PubMed]
Kitagawa, G. The two-filter formula for smoothing and an implementation of the Gaussian-sum smoother. Ann. Inst. Stat. Math. 1994, 46, 605–623. [Google Scholar] [CrossRef]
Kitagawa, G.; Gersch, W. Smoothness Priors Analysis of Time Series; Springer Science & Business Media: Berlin/Heidelberg, Germany, 1996; Volume 116. [Google Scholar]
Karlsson, R.; Gustafsson, F. Particle filtering for quantized sensor information. In Proceedings of the 2005 13th European Signal Processing Conference, Antalya, Turkey, 4–8 September 2005; pp. 1–4. [Google Scholar]
DeGroot, M.H. Optimal Statistical Decisions; Wiley Classics Library, Wiley: Hoboken, NJ, USA, 2005. [Google Scholar]
Borovkov, A.A. Probability Theory; Springer: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
Cohen, H. Numerical Approximation Methods; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
Arasaratnam, I.; Haykin, S.; Elliott, R.J. Discrete-time nonlinear filtering algorithms using Gauss–Hermite quadrature. Proc. IEEE 2007, 95, 953–977. [Google Scholar] [CrossRef]
Arasaratnam, I.; Haykin, S. Cubature kalman smoothers. Automatica 2011, 47, 2245–2250. [Google Scholar] [CrossRef]
Runnalls, A.R. Kullback-Leibler approach to Gaussian mixture reduction. IEEE Trans. Aerosp. Electron. Syst. 2007, 43, 989–999. [Google Scholar] [CrossRef]
Balenzuela, M.P.; Dahlin, J.; Bartlett, N.; Wills, A.G.; Renton, C.; Ninness, B. Accurate Gaussian Mixture Model Smoothing using a Two-Filter Approach. In Proceedings of the 2018 IEEE Conference on Decision and Control (CDC), Miami, FL, USA, 17–19 December 2018; pp. 694–699. [Google Scholar]
Reddy, H.P.; Narasimhan, S.; Bhallamudi, S.M. Simulation and state estimation of transient flow in gas pipeline networks using a transfer function model. Ind. Eng. Chem. Res. 2006, 45, 3853–3863. [Google Scholar] [CrossRef]
Chang, S.I.; Wang, S.J.; Lin, M.S. The Role of Simulation in Control System Design/Modification. In Nuclear Simulation; Springer: Berlin/Heidelberg, Germany, 1987; pp. 205–222. [Google Scholar]
Anderson, B.D.O.; Moore, J.B. Optimal Control: Linear Quadratic Methods; Courier Corporation: Massachusetts, MA, USA, 2007. [Google Scholar]
Gómez, J.C.; Sad, G.D. A State Observer from Multilevel Quantized Outputs. In Proceedings of the 2020 Argentine Conference on Automatic Control (AADECA), Buenos Aires, Argentina, 28–30 October 2020; pp. 1–6. [Google Scholar]
Goodwin, G.C.; Graebe, S.F.; Salgado, M.E. Control System Design; Prentice Hall: Hoboken, NJ, USA, 2001; Volume 240. [Google Scholar]

Figure 1. Diagram of a real dynamic system with quantized data.

Figure 2. State-space model with quantized output.

Figure 3. Representation of the (left) uniform infinite- and (right) finite-level quantizers defined in terms of the quantized values

β_{i}

and the intervals

J_{i}

Figure 3. Representation of the (left) uniform infinite- and (right) finite-level quantizers defined in terms of the quantized values

β_{i}

and the intervals

J_{i}

Figure 4. The shaded area represents the probability of

y_{t}

taking a value

β_{i}

that is equal to the probability of

z_{t}

belonging to set

\in J_{i}

. Here,

P \{\cdot\}

denotes probability.

Figure 4. The shaded area represents the probability of

y_{t}

taking a value

β_{i}

that is equal to the probability of

z_{t}

belonging to set

\in J_{i}

. Here,

P \{\cdot\}

denotes probability.

Figure 5. Example 1. Filtering PDFs at some instants of time for

Δ_{q} = 7

K = 10

, and 100 particles. GT stands for ground truth. QKF, KF, PF, and GSF stand for quantized Kalman filter, Kalman filter, particle filter, and Gaussian sum filter, respectively.

Figure 5. Example 1. Filtering PDFs at some instants of time for

Δ_{q} = 7

K = 10

, and 100 particles. GT stands for ground truth. QKF, KF, PF, and GSF stand for quantized Kalman filter, Kalman filter, particle filter, and Gaussian sum filter, respectively.

Figure 6. MSE between the true and estimated state for

Δ_{q} = \{3, 5, 7\}

. KF, QKF, GSF, and PF stand for the Kalman filter, quantized Kalman filter, Gaussian sum filter, and particle filter, respectively.

Figure 6. MSE between the true and estimated state for

Δ_{q} = \{3, 5, 7\}

. KF, QKF, GSF, and PF stand for the Kalman filter, quantized Kalman filter, Gaussian sum filter, and particle filter, respectively.

Figure 7. Smoothing PDF at time

t = 100

, for

K = \{6, 8, 10\}

, where K increases to the right, and for

Δ_{q} = \{3, 5, 7\}

, where

Δ_{q}

increases upwards. GT stands for the ground truth. GSS, KS, QKS, and PS stand for the Gaussian sum smoother, Kalman smoother, quantized Kalman smoother, and particle smoother, respectively.

Figure 7. Smoothing PDF at time

t = 100

, for

K = \{6, 8, 10\}

, where K increases to the right, and for

Δ_{q} = \{3, 5, 7\}

, where

Δ_{q}

increases upwards. GT stands for the ground truth. GSS, KS, QKS, and PS stand for the Gaussian sum smoother, Kalman smoother, quantized Kalman smoother, and particle smoother, respectively.

Figure 8. Two-tank system.

h_{1}

and

h_{2}

denote the liquid level in Tank 1 and Tank 2, respectively. The liquid flows into Tank 1 at a rate

f_{1}

and out of Tank 2 at a rate

f_{2}

. The quantizer has minimum and maximum values

β_{1} = 2

and

β_{L} = 10

Figure 8. Two-tank system.

h_{1}

and

h_{2}

denote the liquid level in Tank 1 and Tank 2, respectively. The liquid flows into Tank 1 at a rate

f_{1}

and out of Tank 2 at a rate

f_{2}

. The quantizer has minimum and maximum values

β_{1} = 2

and

β_{L} = 10

Figure 9. (Left) The quantized and nonquantized measurements; (right) the MSE between the true and estimated filtered state. KF, QKF, GSF, and PF stand for the Kalman filter, quantized Kalman filter, Gaussian sum filter, and particle filter, respectively.

Table 1. Integral limits and parameters of Theorem 1. ILQ and FLQ mean infinite- and finite-level quantizers, respectively.

	FLQ: $y_{t} = β_{1}$	ILQ: $y_{t} = β_{i}$ with $i \in I$ in (7) FLQ: $y_{t} = β_{i}$ with $i = 2, \dots, L - 1$	FLQ: $y_{t} = β_{L}$
$a_{t}$	$- \infty$	$q_{i - 1} - C x_{t} - D u_{t}$	$q_{L - 1} - C x_{t} - D u_{t}$
$b_{t}$	$q_{1} - C x_{t} - D u_{t}$	$q_{i} - C x_{t} - D u_{t}$	∞
$ς_{t}^{τ}$	$2 ω_{τ} / {(1 + ψ_{τ})}^{2}$	$ω_{τ} (q_{i} - q_{i - 1}) / 2$	$2 ω_{τ} / {(1 + ψ_{τ})}^{2}$
$η_{t}^{τ}$	$- (1 - ψ_{τ}) / (1 + ψ_{τ})$	$ψ_{τ} (q_{i} - q_{i - 1}) / 2$	$(1 - ψ_{τ}) / (1 + ψ_{τ})$
$μ_{t}^{τ}$	$- q_{1}$	$- (q_{i} + q_{i - 1}) / 2$	$- q_{L - 1}$

Table 2. Time in seconds required to perform the smoothing algorithm for the scalar system. Par(K) represents the number of particles (labeled in Figure 7) that produce a similar result to the Gaussian sum smoother with K components. KS, QKS, GSS, and PS stand for the Kalman smoother, quantized Kalman smoother, Gaussian sum smoother, and particle smoother, respectively.

	KS	QKS	GSS			PS
$Δ_{q}$	-	-	$K = 6$	$K = 8$	$K = 10$	Par(6)	Par(8)	Par(10)
7	0.0379	0.1494	0.4158	0.4221	0.4735	2.0420	3.3538	4.0462
5	0.0048	0.0044	0.2695	0.3706	0.4588	3.0074	3.5665	4.0192
3	0.0033	0.0043	0.2609	0.3491	0.4618	4.2471	5.5501	6.3408

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Cedeño, A.L.; Albornoz, R.; Carvajal, R.; Godoy, B.I.; Agüero, J.C. A Two-Filter Approach for State Estimation Utilizing Quantized Output Data. Sensors 2021, 21, 7675. https://doi.org/10.3390/s21227675

AMA Style

Cedeño AL, Albornoz R, Carvajal R, Godoy BI, Agüero JC. A Two-Filter Approach for State Estimation Utilizing Quantized Output Data. Sensors. 2021; 21(22):7675. https://doi.org/10.3390/s21227675

Chicago/Turabian Style

Cedeño, Angel L., Ricardo Albornoz, Rodrigo Carvajal, Boris I. Godoy, and Juan C. Agüero. 2021. "A Two-Filter Approach for State Estimation Utilizing Quantized Output Data" Sensors 21, no. 22: 7675. https://doi.org/10.3390/s21227675

APA Style

Cedeño, A. L., Albornoz, R., Carvajal, R., Godoy, B. I., & Agüero, J. C. (2021). A Two-Filter Approach for State Estimation Utilizing Quantized Output Data. Sensors, 21(22), 7675. https://doi.org/10.3390/s21227675

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

A Two-Filter Approach for State Estimation Utilizing Quantized Output Data

Abstract

1. Introduction

Main Contributions

2. Statement of the Problem

2.1. System Model

2.2. Quantizer Model

2.3. Problem Definition

3. Gaussian Sum Filtering and Smoothing for Quantized Data

3.1. Gaussian Mixture Models

3.2. General Bayesian Framework

3.3. Computing an Explicit Model of p ( y t | x t )

3.4. Understanding the Gaussian Sum Filtering and Smoothing Algorithms

3.5. Gaussian Sum Filtering for Quantized Data

3.6. Computing the State Estimator x ^ t | t from a GMM

3.7. Backward Filtering for Quantized Data

3.8. Smoothing Algorithm with Quantized Data

3.9. Computing the Smoothing Joint PDF p ( x t + 1 , x t | y 1 : N )

4. Numerical Example

4.1. Example 1: First-Order System

4.2. Real-World Application: Tank Liquid Level

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Abbreviations

Appendix A. Technical Lemmata

Appendix B. Quantities to Perform the Gaussian Sum Filter

Appendix C. Quantities to Perform the Backward Filter

Appendix D. Quantities to Perform the Gaussian Sum Smoother

Appendix E. Quantities to Compute p(xt+1, xt|y1:N)

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

3.3. Computing an Explicit Model of $p (y_{t} | x_{t})$

3.6. Computing the State Estimator ${\hat{x}}_{t | t}$ from a GMM

3.9. Computing the Smoothing Joint PDF $p (x_{t + 1}, x_{t} | y_{1 : N})$

Appendix E. Quantities to Compute p(x_t+1, x_t|y_1:N)