Gauss process state-space model optimization algorithm with expectation maximization

Abstract

A Gauss process state-space model trained in a laboratory cannot accurately simulate a nonlinear system in a non-laboratory environment. To solve this problem, a novel Gauss process state-space model optimization algorithm is proposed by combining the expectation–maximization algorithm with the Gauss process Rauch–Tung–Striebel smoother algorithm, that is, the EM-GP-RTSS algorithm. First, a theoretical formulation of the Gauss process state-space model is proposed, which is not found in previous references. Second, a Gauss process state-space model optimization framework with the expectation–maximization algorithm is proposed. In the expectation–maximization algorithm, the unknown system state is considered as the lost data, and the maximization of measurement likelihood function is transformed into that of a conditional expectation function. Then, the Gauss process–assumed density filter algorithm and the Gauss process Rauch–Tung–Striebel smoother algorithm are proposed with the Gauss process state-space model defined in this article, in order to calculate the smoothed distribution in the conditional expectation function. Finally, the Monte Carlo numerical integral method is used to obtain the approximate expression of the conditional expectation function. The simulation results demonstrate that the Gauss process state-space model optimized by the EM-GP-RTSS can simulate the system in the non-laboratory environment better than the Gauss process state-space model trained in the laboratory, and can reach or exceed the estimation accuracy of the traditional state-space model.

Introduction

In many engineering problems, there are a large number of complex nonlinear dynamic systems needed to estimate state with a time series of the sensor measurement. A nonlinear dynamic model should be established according to a nonlinear system, such as a nonlinear state-space model (SSM), consisting of a parametric state transition equation and measurement equation. The Gauss process (GP) model can be used to simulate the nonlinear state transition equation and the measurement equation separately. A non-parametric Gauss process state-space model (GPSSM) can be constructed, if we lack the knowledge of a nonlinear system to build a principled parametric model. For example, GPSSMs are used for target monitoring and location problems in wireless sensor networks.^1,2

In order to train the GPSSM of a nonlinear system, its true system states are required. However, paradoxically, its states are usually unknown and needed to be estimated. At present, some methods and models are proposed to try to overcome training and estimating problems of the GPSSM.

Gauss process latent variable models (GPLVMs) ³ are first proposed to learn a system model with low-dimensional state-space and high-dimensional observation space. Later, GPLVMs are extended to the dynamic system to obtain GPBF-LEARN, which is a framework for learning GP-BayesFilters from weakly labeled training data only. ⁴ With specially tailored Particle Markov Chain Monte Carlo (MC) samplers, a fully Bayesian approach to inference and learning (i.e. state estimation and system identification) for GPSSMs is proposed by marginalizing state transition functions. ⁵ Although this method preserves a non-parametric representation of a system model, it requires that its measurement model should be known. Variational GPSSM is proposed by combining variational Bayes with sequential MC. ⁶ Through imposing a structured Gaussian variational posterior distribution over latent states, an inference mechanism for the GPSSM is proposed in order to address the nonlinear system identification. ⁷ Although the above works have been devoted to the GPSSM training based on the unknown true system states, unfortunately the problem has not been solved theoretically.

Thus, the best way to obtain an accurate GPSSM is still to have to acquire its true system states to train a state transition equation a measurement equation. In order to obtain true system states, it is necessary to build a state detection system with multiple types of sensors, which is expensive and not possible in some cases. When a homogeneous nonlinear system can be placed in a laboratory, such a state detection system can be constructed in the laboratory. After true states are obtained in the laboratory, a GPSSM of the homogeneous nonlinear system can be established.

However, our purpose is to simulate a nonlinear dynamic system in a practical application with GPSSMs, that is, the nonlinear dynamic system in a non-laboratory environment. Because noise and parameters in a non-laboratory environment differ from those of the homogeneous nonlinear system in a laboratory, simulating a homogeneous nonlinear system in a non-laboratory environment with a GPSSM trained in a laboratory could introduce a model bias, which will lead to an overfitting problem of state estimation.

In order to improve the accuracy of a GPSSM in a non-laboratory environment, a reasonable method is to eliminate its model bias by adaptively correcting the hyperparameters of the GPSSM that is trained in a laboratory.

Recently, parameter identification methods, combining the expectation–maximization (EM) algorithm ⁸ with smoothing algorithms, ⁹ are gradually recognized for parametric nonlinear SSMs. Given the input measurements and the output measurements, the EM estimates the parameters in a multivariable bilinear SSM as the robust and gradient-free computation of the maximum likelihood (ML) solution. ¹⁰ A parameter estimation method combining particle method with the EM is proposed to estimate the general SSM of nonlinear systems. ¹¹ Recognizing the unknown noise covariance of a nonlinear system is studied through the EM and the moving horizon estimation. ¹² The EM using the classical extended and ensemble versions of the Kalman smoother is proposed to estimate the model error covariance. ¹³ For nonlinear SSMs with multiplicative unknown parameters, an adaptive dispersion filter algorithm based on the EM is proposed. ¹⁴

Given training data, a GPSSM can be regarded as a special parametric SSM. Based on the maximum likelihood estimator/estimate (MLE) criterion, this article proposes a novel GPSSM optimization algorithm called as the EM-GP-RTSS, combining the EM algorithm with the Gauss process Rauch–Tung–Striebel smoother (GP-RTSS) algorithm, in order to improve the accuracy of GPSSMs simulating the nonlinear dynamic system in a non-laboratory environment.

First, a theoretical formulation of GPSSMs is proposed by expanding one-dimensional output prediction to a multidimensional output prediction in this article. This is an important work, but not seen in previous references. Second, the MLE criterion is used to describe the GPSSM optimization problem as the maximization of the measurement likelihood function. However, the analytical expression of the measurement likelihood function is difficult, and thus the EM is used to replace the difficult-to-maximize measurement likelihood function with the conditional expectation function (CEF). Third, a GPSSM optimization framework based on the EM is proposed in this article. The GP-RTSS algorithm is given to obtain the smoothed estimation distribution of the system state, which is required by the calculation of the CEF. Fourth, in order to calculate the integral of the CEF in the closed form, the MC numerical integral method is used to express an approximate CEF. Fifth, the variable metric method ¹⁴ is adopted for the optimization calculation of the approximate CEF. In the simulation experiment, a noise pendulum model is used as a representative of the nonlinear dynamic system to verify the effectiveness of the EM-GP-RTSS algorithm proposed in this article.

This article is organized as follows. In section “Background on GP,” we briefly review the background on GPs. In section “Optimization algorithm for GPSSM,” we describe the development of the GPSSM optimization algorithm in detail. In section “Simulation experiment,” the experimental results based on the noise pendulum model are presented to demonstrate the performance of the EM-GP-RTSS. Finally, the conclusion is drawn in section “Conclusion and future works.”

Background on GP

In this section, a standard GP model and a noise input prediction model are given. This section is the mathematical foundation for the development of the GPSSM optimization algorithm.

Standard GP model

A GP is a non-parametric tool for learning regression functions from the sample data. ¹⁵ Specifically, a GP represents a random process, whose outputs of sampling points obey to a joint Gaussian distribution. A GP can be fully described by its mean function and covariance function. Suppose there is a training data set from the following noise process

$${\beta }_i=\psi \left({\mathit{\alpha }}_i\right)+\epsilon$$

(1)

where ${\mathit{\alpha }}_i$ is an n-dimensional input column vector, ${\beta }_i$ is a scalar output, $\epsilon$ is a random variable that obeys a Gaussian white noise $N(\epsilon |0,{\sigma }^2)$ , where $\sigma$ is its standard deviation. Through denoting the matrix $\mathit{A}=[{\mathit{\alpha }}_1,{\mathit{\alpha }}_2,\dots ,{\mathit{\alpha }}_N]$ and the column vector $\mathit{\beta }=({\beta }_1,{\beta }_2,\dots ,{\beta }_N{)}^T$ , the data set D can be restated as $D=〈\mathit{A},\mathit{\beta }〉$ . The joint distribution of the noise output $\mathit{\beta }$ is a function of $\mathit{A}$ . This distribution is assumed as a zero-mean multidimensional Gaussian distribution

$$p\left(\mathit{\beta }\right)=N(\mathit{0},\mathit{K}(\mathit{A},\mathit{A})+{\sigma }^2\mathit{I})$$

(2)

where $\mathit{K}(\mathit{A},\mathit{A})$ is the kernel matrix, and its element of the $i\mathrm{th}$ row and the $j\mathrm{th}$ column is denoted as ${\mathit{K}}_{\mathrm{ij}}=k({\mathit{\alpha }}_i,{\mathit{\alpha }}_j)$ . The kernel function of GPSSMs usually selects the squared exponential (SE) function

$$k({\mathit{\alpha }}_i,{\mathit{\alpha }}_j)={\sigma }_{\psi }^2\exp (-\frac{1}{2}{({\mathit{\alpha }}_i-{\mathit{\alpha }}_j)}^T{\mathit{W}}_{\psi }^{-1}({\mathit{\alpha }}_i-{\mathit{\alpha }}_j))$$

(3)

where ${\sigma }_{\psi }^2$ is a signal variance. The diagonal matrix ${\mathit{W}}_{\psi }$ consists of the length scales for each dimension output. ${\sigma }_{\psi }$ , $\sigma$ , and ${\mathit{W}}_{\psi }$ are the hyperparameters of the GP. Given a training data set $〈\mathit{A},\mathit{\beta }〉$ and a noiseless test input ${\mathit{\alpha }}_{*}$ , the Gaussian distribution of the output ${\beta }_{*}$ can be obtained by the standard GP prediction model, whose mean is

$$G{P}_{\rho }({\mathit{\alpha }}_{*},〈\mathit{A},\mathit{\beta }〉)={\mathit{k}}_{*}^T{[\mathit{K}(\mathit{A},\mathit{A})+{\sigma }^2\mathit{I}]}^{-1}\mathit{\beta }$$

(4)

and whose variance is

$$G{P}_{\Sigma }({\mathit{\alpha }}_{*},〈\mathit{A},\mathit{\beta }〉)=k({\mathit{\alpha }}_{*},{\mathit{\alpha }}_{*})-{\mathit{k}}_{*}^T{[\mathit{K}(\mathit{A},\mathit{A})+{\sigma }^2\mathit{I}]}^{-1}{\mathit{k}}_{*}$$

(5)

where ${\mathit{k}}_{*}$ is the kernel vector between ${\mathit{\alpha }}_{*}$ and $\mathit{A}$ , denoted as $(k({\mathit{\alpha }}_{*},{\mathit{\alpha }}_1),k({\mathit{\alpha }}_{*},{\mathit{\alpha }}_2),\dots ,k({\mathit{\alpha }}_{*},{\mathit{\alpha }}_N){)}^T$ .

Noise input prediction model

If the test input ${\mathit{\alpha }}_{*}$ is a random Gaussian noise, the prediction distribution of corresponding output ${\beta }_{*}$ cannot be directly given by formulas (4) and (5). Assume that ${\mathit{\alpha }}_{*}$ follows the Gaussian distribution $N({\mathit{\alpha }}_{*}|{\hat{\mathit{\alpha }}}_{*},{\mathit{P}}_{*})$ , where ${\hat{\mathit{\alpha }}}_{*}$ is a mean vector and ${\mathit{P}}_{*}$ is a covariance matrix. The prediction problem corresponds to calculating the probability distribution

$$p\left({\beta }_{*}\right)=\int \int p({\beta }_{*}|{\mathit{\alpha }}_{*},\epsilon )p\left({\mathit{\alpha }}_{*}\right)p\left(\epsilon \right)d{\mathit{\alpha }}_{*}d\epsilon$$

(6)

The mean and variance of $p({\beta }_{*}|{\mathit{\alpha }}_{*},\epsilon )$ are calculated by formulas (4) and (5), respectively. Fortunately, when the SE kernel function is used, the mean ${\hat{\beta }}_{*}$ and the variance ${\sigma }_{\beta }^2$ of $p\left({\beta }_{*}\right)$ can be calculated analytically

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\beta }}_{*}={\mathrm{E}}_{\epsilon }\left[{\mathrm{E}}_{{\mathit{\alpha }}_{*}}\left[{\mathrm{E}}_{\psi }[\psi \left({\mathit{\alpha }}_{*}\right)+\epsilon |{\mathit{\alpha }}_{*},\epsilon ]\right]\right]}\\ \multicolumn{1}{c}{=\int \int G{P}_{\mu }({\mathit{\alpha }}_{*},〈\mathit{A},\mathit{\beta }〉)N\left({\mathit{\alpha }}_{*}\right|{\hat{\mathit{\alpha }}}_{*},{\mathit{P}}_{*})N(\epsilon |0,{\sigma }^2)d{\mathit{\alpha }}_{*}d\epsilon }\\ \multicolumn{1}{c}{={\mathbf{\lambda }}^T\mathbf{\gamma }}\end{array}$$

(7)

where $\mathbf{\lambda }={[\mathit{K}\left(\mathit{A},\mathit{A}\right)+{\sigma }^2\mathit{I}]}^{-1}\mathit{\beta }$ and $\mathbf{\gamma }=({\gamma }_1,{\gamma }_2,\dots ,{\gamma }_N)$ . $\mathrm{E}[·]$ and $\mathrm{E}[·|·]$ represent the expectation operation and the conditional expectation operation, respectively.

Let $i=1,2,\dots ,N$ , and ${\gamma }_i$ is

$$\begin{array}{c}\multicolumn{1}{c}{{\gamma }_i=\int k({\mathit{\alpha }}_i,{\mathit{\alpha }}_{*})N\left({\mathit{\alpha }}_{*}\right|{\hat{\mathit{\alpha }}}_{*},{\mathit{P}}_{*})d{\mathit{\alpha }}_{*}}\\ \multicolumn{1}{c}{={\sigma }_{\psi }^2{|{\mathit{P}}_{*}{\mathit{W}}_{\psi }^{-1}+\mathit{I}|}^{-\frac{1}{2}}}\\ \multicolumn{1}{c}{\exp (-\frac{1}{2}{({\mathit{\alpha }}_i-{\hat{\mathit{\alpha }}}_{*})}^T{({\mathit{W}}_{\psi }+{\mathit{P}}_{*})}^{-1}({\mathit{\alpha }}_i-{\hat{\mathit{\alpha }}}_{*}))}\end{array}$$

(8)

The variable ${\gamma }_i$ is the expectation of $k({\mathit{\alpha }}_i,{\mathit{\alpha }}_{*})$ on ${\mathit{\alpha }}_{*}$ . ${\sigma }_{\beta }^2$ is

$$\begin{array}{c}{\sigma }_{\beta }^2=\int {\int \int ({\beta }_{*}-{\hat{\beta }}_{*})}^{2}p({\beta }_{*}|{\mathit{\alpha }}_{*},\epsilon )p\left({\mathit{\alpha }}_{*}\right)p\left(\epsilon \right)d{\beta }_{*}d{\mathit{\alpha }}_{*}d\epsilon \hfill \\ =\int [\int {({\beta }_{*}-{\hat{\beta }}_{*})}^{2}p({\beta }_{*}|{\mathit{\alpha }}_{*},\epsilon )d{\beta }_{*}]p\left({\mathit{\alpha }}_{*}\right)p\left(\epsilon \right)d{\mathit{\alpha }}_{*}d\epsilon \hfill \\ =\int [G{P}_{\mu }{({\mathit{\alpha }}_{*},〈\mathit{A},\mathit{\beta }〉)}^2+G{P}_{\Sigma }({\mathit{\alpha }}_{*},〈\mathit{A},\mathit{\beta }〉)+{\sigma }^2-{\hat{\beta }}_{*}^2]\hfill \\ \multicolumn{1}{c}{p\left({\mathit{\alpha }}_{*}\right)d{\mathit{\alpha }}_{*}}\\ ={\mathbf{\lambda }}^T\mathit{\Phi }\mathbf{\lambda }+{\sigma }_{\psi }^2-\mathrm{tr}\left({(\mathit{K}(\mathit{A},\mathit{A})+{\sigma }^2\mathit{I})}^{-1}\mathit{\Phi }\right)-{\hat{\beta }}_{*}^2+{\sigma }^2\hfill \end{array}$$

(9)

where $\mathrm{tr}(·)$ is the trace operation. The element of the $i\mathrm{th}$ row and the $j\mathrm{th}$ column in $\mathit{\Phi }$ is

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{\Phi }}_{\mathrm{ij}}=\frac{k({\mathit{\alpha }}_i,{\hat{\mathit{\alpha }}}_{*})k({\mathit{\alpha }}_j,{\hat{\mathit{\alpha }}}_{*})}{{|2{\mathit{P}}_{*}{\mathit{W}}_{\psi }^{-1}+\mathit{I}|}^{-\frac{1}{2}}}}\\ \multicolumn{1}{c}{\times \exp \left({({\tilde{\mathit{\alpha }}}_{\mathrm{ij}}-{\hat{\mathit{\alpha }}}_{*})}^T{(\frac{1}{2}{\mathit{W}}_{\psi }+{\mathit{P}}_{*})}^{-1}{\mathit{P}}_{*}{\mathit{W}}_{\psi }^{-1}({\tilde{\mathit{\alpha }}}_{\mathrm{ij}}-{\hat{\mathit{\alpha }}}_{*})\right)}\end{array}$$

(10)

where ${\tilde{\mathit{\alpha }}}_{\mathrm{ij}}=1/2({\mathit{\alpha }}_i+{\mathit{\alpha }}_j)$ . Thus, $p({\beta }_{*})$ can be approximated by $N({\beta }_{*}|{\hat{\beta }}_{*},{\sigma }_{\beta }^2)$ in the closed form.

Optimization algorithm for GPSSM

In this section, the GPSSM optimization algorithm is developed with the EM and the GP-RTSS. The whole framework flowchart of the GPSSM optimization algorithm is shown in Figure 1. First, the training data are obtained to train an original GPSSM to simulate a nonlinear system in a laboratory. Second, the measurements of nonlinear dynamic systems are obtained in a non-laboratory environment. The EM is used to optimize the GPSSM of homogeneous system in the laboratory. Since calculating the CEF in the Expectation step (E-step) requires the smoothed distribution of system states, the GP-RTSS algorithm is introduced into the EM optimization framework. In the Maximization step (M-step), the CEF is maximized to obtain the iterative hyperparameters, which means to adjust the GPSSM at the same time. The E-step and the M-step are repeated alternately until convergence occurs, and the optimized GPSSM and the estimated states will be acquired at the end.

Figure 1.

The flowchart of the GPSSM optimization.

GPSSM

This section gives a theoretical formulation of GPSSMs. The concept of GPSSMs is already proposed in Ko and Fox,^4,16 and Deisenroth et al., ¹⁷ but a theoretical formulation of GPSSMs has not been found. Consider the traditional SSM of a discrete-time nonlinear system

$$\begin{array}{c}{\mathit{x}}_{k+1}=\mathit{f}\left({\mathit{x}}_k\right)+{\mathit{w}}_k\hfill \\ {\mathit{y}}_{k+1}=\mathit{h}\left({\mathit{x}}_{k+1}\right)+{\mathit{v}}_{k+1}\hfill \end{array}$$

(11)

where $k\in \{0,1,2,\dots \}$ is an index of the time step, ${\mathit{x}}_k\in {\mathbb{R}}^n$ is a state vector, ${\mathit{y}}_k\in {\mathbb{R}}^m$ is a measurement vector, ${\mathit{w}}_k$ and ${\mathit{v}}_k$ are a zero-mean Gaussian white noise, $n,m\in \mathbb{N}$ . The functions $\mathit{f}$ and $\mathit{h}$ are a state transition function and a measurement function, respectively. If the dynamic system is well understood, a traditional modeling method is to construct the parametric expressions of $\mathit{f}$ and $\mathit{h}$ . If not, it is very difficult to establish an appropriate parametric model. However, if the GP prior distributions of $\mathit{f}$ and $\mathit{h}$ are set separately, the GPSSM can be used to simulate the nonlinear SSM described in equation (11).

The operating mechanism of the dynamic system is assumed to be constant, that is the corresponding SSM of the system is unchanged in a time window $[t_0,t_0+\mathrm{\Delta }t\times L]$ , where $t_0$ is the left edge of the time window, $t_0+\mathrm{\Delta }t\times L$ is the right edge of the time window, L is a sample number, and $\mathrm{\Delta }t$ is a sampling period. For ease of writing, set $t_0=0\mathrm{s}$ and $\mathrm{\Delta }t=1\mathrm{s}$ in this article.

The matrix ${\mathit{X}}_{0:L}=[{\mathit{x}}_0,{\mathit{x}}_1,\dots ,{\mathit{x}}_L]$ is a sequence of states in the state-space, and the matrix ${\mathit{Y}}_{1:L}=[{\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_L]$ is a sequence of measurements for ${\mathit{X}}_{1:L}$ . Because the output of GP is a scalar, $\mathit{f}$ and $\mathit{h}$ need to be decomposed into several one-dimensional functions, namely, $\mathit{f}=(f_1,f_2,\dots ,f_n{)}^T$ and $\mathit{h}=(h_1,h_2,\dots ,h_m{)}^T$ . The SSM (equation (11)) can be expressed as follows

$$\{\begin{array}{c}{\mathit{x}}_{k+1}^{(1)}=f_1\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(1)}\hfill \\ \multicolumn{1}{c}{{\mathit{x}}_{k+1}^{(2)}=f_2\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(2)}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{x}}_{k+1}^{(n)}=f_n\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(n)}}\end{array};\{\begin{array}{c}{\mathit{y}}_{k+1}^{(1)}=h_1\left({\mathit{x}}_{k+1}\right)+{\mathit{v}}_{k+1}^{(1)}\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_{k+1}^{(2)}=h_2\left({\mathit{x}}_{k+1}\right)+{\mathit{v}}_{k+1}^{(2)}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{y}}_{k+1}^{(m)}=h_m\left({\mathit{x}}_{k+1}\right)+{\mathit{v}}_{k+1}^{(n)}}\end{array}$$

(12)

The set $D_{f_i}=〈{\mathit{X}}_{0:L-1},{\mathit{x}}_{1:L}^{(i)}〉$ is denoted as the training data set of $f_i$ , where ${\mathit{X}}_{0:L-1}$ is the training input data, ${\mathit{x}}_{1:L}^{(i)}=({\mathit{x}}_1^{(i)},{\mathit{x}}_2^{(i)},\dots ,{\mathit{x}}_L^{(i)}{)}^T$ is the training output data, ${\mathit{x}}^{(i)}$ is the $i\mathrm{th}$ dimension element of $\mathit{x}$ . Similarly, $D_{h_i}=〈{\mathit{X}}_{1:L},{\mathit{y}}_{1:L}^{(i)}〉$ is denoted as the training data set of $h_i$ , where ${\mathit{y}}_{1:L}^{(i)}=({\mathit{y}}_1^{(i)},{\mathit{y}}_2^{(i)},\dots ,{\mathit{y}}_L^{(i)}{)}^T$ . The prior distribution of the output training data ${\mathit{x}}_{1:L}^{(i)}$ and ${\mathit{y}}_{1:L}^{(i)}$ can be denoted as two zero-mean multidimensional Gaussian distributions, respectively

$$p\left({\mathit{x}}_{1:L}^{(i)}\right)=N(0,{\mathit{K}}_{f_i}({\mathit{X}}_{0:L-1},{\mathit{X}}_{0:L-1})+{\sigma }_{\mathrm{xi}}^2\mathit{I})$$

(13)

$$p\left({\mathit{y}}_{1:L}^{(i)}\right)=N(0,{\mathit{K}}_{h_i}({\mathit{X}}_{1:L},{\mathit{X}}_{1:L})+{\sigma }_{\mathrm{yi}}^2\mathit{I})$$

(14)

where ${\mathit{K}}_{f_i}({\mathit{X}}_{0:L-1},{\mathit{X}}_{0:L-1})$ and ${\mathit{K}}_{h_i}({\mathit{X}}_{1:L},{\mathit{X}}_{1:L})$ are the kernel matrixes abbreviated as ${\mathit{K}}_{f_i}$ and ${\mathit{K}}_{h_i}$ . The SE kernel functions of ${\mathit{K}}_{f_i}$ and ${\mathit{K}}_{h_i}$ are denoted as, respectively

$$k_{f_i}(\mathit{x},{\mathit{x}}^{\prime })={\sigma }_{f_i}^2\exp (-\frac{1}{2}{(\mathit{x}-{\mathit{x}}^{\prime })}^T{\mathit{W}}_{f_i}^{-1}(\mathit{x}-{\mathit{x}}^{\prime }))$$

(15)

$$k_{h_i}(\mathit{x},{\mathit{x}}^{\prime })={\sigma }_{h_i}^2\exp (-\frac{1}{2}{(\mathit{x}-{\mathit{x}}^{\prime })}^T{\mathit{W}}_{h_i}^{-1}(\mathit{x}-{\mathit{x}}^{\prime }))$$

(16)

The above formulas (12)–(16) define a formulated GPSSM, which can simulate the traditional nonlinear SSM described in formula (11).

Because the system state is unknown, the state ${\mathit{x}}_k$ is assumed to obey a prior Gaussian distribution $N({\mathit{x}}_k|{\hat{\mathit{x}}}_k,{\mathit{P}}_k)$ , where ${\hat{\mathit{x}}}_k$ is the mean vector and ${\mathit{P}}_k$ is the covariance matrix. The state ${\mathit{x}}_k$ is the noise test input of the GPSSM, and then ${\mathit{x}}_{k+1}$ and ${\mathit{y}}_k$ are denoted as the state transition output and the measurement output of the GPSSM, respectively. The prediction distribution of ${\mathit{x}}_{k+1}$ and ${\mathit{y}}_k$ are assumed as Gaussian distributions $N({\mathit{x}}_{k+1}|{\overline{\mathit{x}}}_{k+1},{\mathit{Q}}_{k+1})$ and $N({\mathit{y}}_k|{\overline{\mathit{y}}}_k,{\mathit{R}}_k)$ , respectively, where ${\overline{\mathit{x}}}_{k+1}$ and ${\overline{\mathit{y}}}_k$ are mean vectors, and ${\mathit{Q}}_{k+1}$ and ${\mathit{R}}_k$ are covariance matrixes. Each elements in ${\overline{\mathit{x}}}_{k+1}$ and ${\overline{\mathit{y}}}_k$ can be calculated by formulas (6)–(8) in the closed form, and ${\mathit{Q}}_{k+1}^{(i,j)}(i,j=1,2,\dots ,n)$ , which is the element of $i\mathrm{th}$ row and the jth column in ${\mathit{Q}}_{k+1}$ , can be obtained by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{Q}}_{k+1}^{(i,j)}=\int ({\mathit{x}}_{k+1}^{(i)}-{\overline{\mathit{x}}}_{k+1}^{(i)})({\mathit{x}}_{k+1}^{(j)}-{\overline{\mathit{x}}}_{k+1}^{(j)})}\\ \multicolumn{1}{c}{p({\mathit{x}}_{k+1}^{(i)}|{\mathit{x}}_k,{\mathit{w}}_k^{(i)})}\\ \multicolumn{1}{c}{·p({\mathit{x}}_{k+1}^{(j)}|{\mathit{x}}_k,{\mathit{w}}_k^{(j)})p\left({\mathit{x}}_k\right)p\left({\mathit{w}}_k^{(i)}\right)p\left({\mathit{w}}_k^{(j)}\right)}\\ \multicolumn{1}{c}{d{\mathit{x}}_{k+1}^{(i)}d{\mathit{x}}_{k+1}^{(j)}d{\mathit{x}}_{k}d{\mathit{w}}_k^{(i)}d{\mathit{w}}_k^{(j)}}\end{array}$$

(17)

If $i=j$ , then ${\mathit{Q}}_{k+1}^{(i,j)}$ is the diagonal element of ${\mathit{Q}}_{k+1}$ , namely, the variance of ${\mathit{x}}_{k+1}^{(i)}$ calculated by formula (9). If $i\ne j$ , it is calculated by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{Q}}_{k+1}^{(i,j)}=\int (f_i\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(i)}-{\overline{\mathit{x}}}_{k+1}^{(i)})(f_j\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(j)}-{\overline{\mathit{x}}}_{k+1}^{(j)})}\\ \multicolumn{1}{c}{p(f_i\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(i)}|{\mathit{x}}_k,{\mathit{w}}_k^{(i)})}\\ \multicolumn{1}{c}{·p(f_j\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(j)}|{\mathit{x}}_k,{\mathit{w}}_k^{(j)})p\left({\mathit{x}}_k\right)p\left({\mathit{w}}_k^{(i)}\right)p\left({\mathit{w}}_k^{(j)}\right)d{\mathit{x}}_{k+1}^{(i)}d{\mathit{x}}_{k+1}^{(j)}}\\ \multicolumn{1}{c}{d{\mathit{x}}_{k}d{\mathit{w}}_k^{(i)}d{\mathit{w}}_k^{(j)}}\\ \multicolumn{1}{c}{=\int (f_i\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(i)}-{\overline{\mathit{x}}}_{k+1}^{(i)})(f_j\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(j)}-{\overline{\mathit{x}}}_{k+1}^{(j)})}\\ \multicolumn{1}{c}{p(f_i\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(i)}|{\mathit{x}}_k,{\mathit{w}}_k^{(i)})}\\ \multicolumn{1}{c}{·p(f_j\left({\mathit{x}}_k\right)+{\mathit{w}}_k^{(j)}|{\mathit{x}}_k,{\mathit{w}}_k^{(j)})p\left({\mathit{x}}_k\right)p\left({\mathit{w}}_k^{(i)}\right)p\left({\mathit{w}}_k^{(j)}\right)d{\mathit{x}}_{k+1}^{(i)}d{\mathit{x}}_{k+1}^{(j)}}\\ \multicolumn{1}{c}{d{\mathit{x}}_{k}d{\mathit{w}}_k^{(i)}d{\mathit{w}}_k^{(j)}}\\ \multicolumn{1}{c}{=\int G{P}_{\rho }({\mathit{x}}_k,D_{f_i})G{P}_{\rho }({\mathit{x}}_k,D_{f_j})p\left({\mathit{x}}_k\right)d{\mathit{x}}_k-{\overline{\mathit{x}}}_{k+1}^{(i)}{\overline{\mathit{x}}}_{k+1}^{(j)}}\\ \multicolumn{1}{c}{=\int k_{f_i}({\mathit{x}}_k,{\mathit{X}}_{0:L-1}){\mathbf{\lambda }}_{f_i}{k}_{f_j}({\mathit{x}}_k,{\mathit{X}}_{0:L-1}){\mathbf{\lambda }}_{f_j}p\left({\mathit{x}}_k\right)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{={\mathbf{\lambda }}_{f_i}^T\int k_{f_i}({\mathit{x}}_k,{\mathit{X}}_{0:L-1})k_{f_j}({\mathit{x}}_k,{\mathit{X}}_{0:L-1})p\left({\mathit{x}}_k\right)d{\mathit{x}}_k{\mathbf{\lambda }}_{f_j}}\\ \multicolumn{1}{c}{={\mathbf{\lambda }}_{f_i}^T{\tilde{\Phi }\mathbf{\lambda }}_{f_j}}\end{array}$$

(18)

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\Phi }}^{(p,q)}={\sigma }_{f_i}^2{\sigma }_{f_j}^2{|({\mathit{W}}_{f_i}^{-1}+{\mathit{W}}_{f_j}^{-1}){\mathit{P}}_k+\mathit{I}|}^{-\frac{1}{2}}}\\ \multicolumn{1}{c}{·\exp (-\frac{1}{2}{({\mathit{x}}_p-{\mathit{x}}_q)}^T{({\mathit{W}}_{f_i}+{\mathit{W}}_{f_j})}^{-1}({\mathit{x}}_p-{\mathit{x}}_q))}\\ \multicolumn{1}{c}{\exp (-\frac{1}{2}{({\mathit{z}}_{\mathrm{pq}}-{\hat{\mathit{x}}}_k)}^T{\mathbf{\Lambda }}^{-1}({\mathit{z}}_{\mathrm{pq}}-{\hat{\mathit{x}}}_k))}\end{array}$$

(19)

where ${\mathbf{\lambda }}_{f_i}={({\mathit{K}}_{f_i}+{\sigma }_{xi}^2\mathit{I})}^{-1}{\mathit{x}}_{1:L}^{(i)}$ , $\mathit{\Lambda }={({\mathit{W}}_{f_i}^{-1}+{\mathit{W}}_{f_j}^{-1})}^{-1}+{\mathit{P}}_k$ , ${\mathit{z}}_{\mathrm{pq}}={\mathit{W}}_{f_j}({\mathit{W}}_{f_i}+{\mathit{W}}_{f_j}{)}^{-1}{\mathit{x}}_p+{\mathit{W}}_{f_i}({\mathit{W}}_{f_i}+{\mathit{W}}_{f_j}{)}^{-1}{\mathit{x}}_q$ , and $p,q=0,1,\dots ,L-1$ .Thus, the mean vector and covariance matrix of the predicted distribution of ${\mathit{x}}_{k+1}$ and ${\mathit{y}}_k$ can be calculated accurately. In addition, the condition prediction distributions of ${\mathit{x}}_{k+1}$ and ${\mathit{y}}_k$ given a noise-free test input ${\mathit{x}}_k$ are calculated, respectively, by

$$\begin{array}{c}\multicolumn{1}{c}{p\left({\mathit{x}}_{k+1}|{\mathit{x}}_k\right)=\underset{i=1}{\overset{n}{\Pi }}}\\ \multicolumn{1}{c}{N({\mathit{x}}_{k+1}^{(i)}|G{P}_{\rho }({\mathit{x}}_k,D_{f_i}),G{P}_{\Sigma }({\mathit{x}}_k,D_{f_i})+{\sigma }_{\mathrm{xi}}^2)}\end{array}$$

(20)

$$p\left({\mathit{y}}_k|{\mathit{x}}_k\right)=\underset{i=1}{\overset{m}{\Pi }}N({\mathit{y}}_k^{(i)}|G{P}_{\rho }({\mathit{x}}_k,D_{h_i}),G{P}_{\Sigma }({\mathit{x}}_k,D_{h_i})+{\sigma }_{\mathrm{yi}}^2)$$

(21)

Optimization algorithm framework

A GPSSM is theoretically expressed in the previous sections. In order to establish an instantiated GPSSM for a specific nonlinear system, the training data set of the state transition equation and the measurement equation should be given, that is, true states ${\mathit{X}}_{0:L}$ and corresponding measurements ${\mathit{Y}}_{1:L}$ need to be known. In a non-laboratory environment, however, ${\mathit{X}}_{0:L}$ is often unknown and needs to be estimated.

An accepted method is to build a state detection system equipped with multiple types of sensors in a laboratory, by which true states can be measured to train a GPSSM. The GPSSM obtained in the laboratory is recorded as the GPSSM* in this article. The GPSSM* can obtain a satisfactory state estimation accuracy in a laboratory, but there is an overfitting problem of state estimation with the GPSSM* in a non-laboratory environment. This is because the environment of non-laboratory and laboratory are different, and then, the homogeneous nonlinear system is affected by several factors. Main differences are system noise, sensor noise, and system parameters. Then, the use of the GPSSM* in a non-laboratory environment would introduce a model bias. The lack of robustness is a drawback of the GPSSM*, but the GPSSM* simulates the basic mechanism of the homogeneous nonlinear system in the non-laboratory environment.

In order to solve this problem, the core idea of this article is to optimize the hyperparameters of the GPSSM* to obtain an approximate optimal model for a nonlinear system in a non-laboratory environment. Its state and measurement sequences obtained in a laboratory are denoted as ${\mathit{X}}_{0:L}^{*}$ and ${\mathit{Y}}_{1:L}^{*}$ , respectively, and the corresponding training data set is denoted as $D^{*}$ . $\{{\sigma }_{f_i},{\mathit{W}}_{f_j},{\sigma }_{\mathrm{xi}}{\}}_{i=1}^n$ and $\{{\sigma }_{h_i},{\mathit{W}}_{h_j},{\sigma }_{\mathrm{yi}}{\}}_{i=1}^m$ of the GPSSM are collectively referred to $\theta$ . Given $D^{*}$ , the MLE of $\theta$ denoted as ${\theta }^{*}$ is acquired by maximizing the measurement likelihood function. ¹⁸ If $D^{*}$ and ${\theta }^{*}$ are given, the GPSSM* is determined and denoted as $G({\theta }^{*},D^{*})$ . A general GPSSM is denoted as $G(\theta ,D)$ . A measurement sequence of a nonlinear system ${\mathit{Y}}_{1:L}$ is obtained in a non-laboratory environment. Based on the MLE criterion, there is an optimal GPSSM $G(\stackrel{⌢}{\theta },\stackrel{⌢}{D})$ that makes the log-likelihood function reach the maximum

$$G(\stackrel{⌢}{\theta },\stackrel{⌢}{D})=\underset{G\left(D,\theta \right)}{\arg \max }\log p\left({\mathit{Y}}_{1:L}|G\left(\theta ,D\right)\right)$$

(22)

Because of the unknown $\hat{D}$ , $G(\stackrel{⌢}{\theta },\stackrel{⌢}{D})$ cannot be obtained. However, an approximate optimal GPSSM $G(\widehat{\theta },D^{\ast })$ can be obtained by correcting ${\theta }^{*}$ of $G({\theta }^{*},D^{*})$

$$\hat{\theta }=\underset{\theta }{\arg max}\log p\left({\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)$$

(23)

Given $D^{*}$ , $G(\theta ,D^{*})$ can be regarded as a special parametric SSM. Formula (23) expresses a hyperparameter optimization problem by the MLE criterion. As the measurement log-likelihood function in formula (23) is difficult to express analytically, the numerical search algorithm cannot be adopted to maximize the function directly.

The unknown system state sequence is denoted as ${\mathit{X}}_{0:L}$ . By the Bayes theorem, the log-likelihood function $\log p({\mathit{Y}}_{1:L}|G(\theta ,D^{*}))$ can be decomposed into

$$\begin{array}{c}\multicolumn{1}{c}{\log p\left({\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)=\log p\left({\mathit{X}}_{0:L},{\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)}\\ \multicolumn{1}{c}{-\log p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G(\theta ,D^{*}))}\end{array}$$

(24)

The hyperparameter $\theta$ at the $\mu \mathrm{th}$ EM iteration is denoted as ${\theta }_{\mu }$ and then perform the expectation operation to formula (24) on $p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))$

$$\begin{array}{c}\multicolumn{1}{c}{\log p\left({\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)=\int p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))}\\ \multicolumn{1}{c}{\log p\left({\mathit{X}}_{0:L},{\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)d{\mathit{X}}_{0:L}}\\ \multicolumn{1}{c}{-\int p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))\log p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G(\theta ,D^{*}))}\\ \multicolumn{1}{c}{d{\mathit{X}}_{0:L}}\end{array}$$

(25)

Let

$$\begin{array}{c}\multicolumn{1}{c}{l(\theta ,{\theta }_{\mu })=\int p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))}\\ \multicolumn{1}{c}{\log p\left({\mathit{X}}_{0:L},{\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)d{\mathit{X}}_{0:L}}\end{array}$$

(26)

Then

$$\begin{array}{c}\log p\left({\mathit{Y}}_{1:L}|G(\theta ,D^{*})\right)-\log p\left({\mathit{Y}}_{1:L}|G({\theta }_{\mu },D^{*})\right)\hfill \\ =l(\theta ,{\theta }_{\mu })-l({\theta }_{\mu },{\theta }_{\mu })+\int p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))\hfill \\ \multicolumn{1}{c}{\log \frac{p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))}{p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G(\theta ,D^{*}))}d{\mathit{X}}_{0:L}}\end{array}$$

(27)

where $l(\theta ,{\theta }_{\mu })$ is a CEF. The last item at the right end of equation (27) is a Kullback–Leibler information distance from $p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G(\theta ,D^{*}))$ to $p({\mathit{X}}_{0:L}|{\mathit{Y}}_{1:L},G({\theta }_{\mu },D^{*}))$ , which proves to be non-negative.

For $\forall {\theta }_{\mu +1}$ , if $l({\theta }_{\mu +1},{\theta }_{\mu })>l({\theta }_{\mu },{\theta }_{\mu })$ , then $p({\mathit{Y}}_{1:L}|G({\theta }_{\mu +1},D^{\ast }))>p({\mathit{Y}}_{1:L}|G({\theta }_{\mu },D^{\ast }))$ . Given $G(\theta ,D^{*})$ and ${\mathit{Y}}_{1:L}$ , an EM-based GPSSM optimization framework is proposed as Framework 1.

• Framework 1 (EM-based GPSSM optimization framework)

• Step 1. Initialize ${\theta }_{\mu }={\theta }^{*}$ and $\mu =0$ ;

• Step 2. Calculate $l(\theta ,{\theta }_{\mu })$ in the E-step;

• Step 3. Maximize ${\mathbf{\theta }}_{\mu +1}=\arg \underset{\theta }{max}l(\mathbf{\theta },{\mathbf{\theta }}_{\mu })$ in the M-step;

• Step 4. Judge the convergence. If not, $\mu =\mu +1$ and return to step 2; otherwise, end.

In the next section, the GPSSM optimization algorithm is developed based on Framework 1.

Expectation step

The first problem in Framework 1 is the specific expression of $l(\theta ,{\theta }_{\mu })$ . For ease of writing, $p(·|G(\theta ,D^{*}))$ is recorded as $p_{\theta }(·)$ . According to the Bayesian conditional probability theorem and the Markov property, $l(\theta ,{\theta }_{\mu })$ can be decomposed as follows

$$l(\theta ,{\theta }_{\mu })=l_0+l_1+l_2$$

(28)

$$l_0=\int p_{{\theta }_{\mu }}({\mathit{x}}_0|{\mathit{Y}}_{1:L})\log p_{\theta }\left({\mathit{x}}_0\right)d{\mathit{x}}_0$$

(29)

$$l_1=\sum _{i=0}^{L-1}\int \int p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})\log p_{\theta }\left({\mathit{x}}_{i+1}|{\mathit{x}}_i\right)d{\mathit{x}}_{i}d{\mathit{x}}_{i+1}$$

(30)

$$l_2=\sum _{i=1}^L\int \int p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})\log p_{\theta }\left({\mathit{y}}_i|{\mathit{x}}_i\right)d{\mathit{x}}_i$$

(31)

where $p_{\theta }({\mathit{x}}_0)$ is a probability distribution of initial state ${\mathit{x}}_0$ , which is assumed to obey the Gaussian distribution with an unknown mean vector ${\hat{\mathit{x}}}_0$ and covariance matrix ${\mathit{P}}_0$ , which are also unknown parameters that are included in $\theta$ .

From formulas (29)–(31), $l(\theta ,{\theta }_{\mu })$ needs to calculate the smoothed distribution $p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})$ and the joint smoothed distribution $p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})$ . When the EM iteratively optimizes $\theta$ , it also estimates the system states. A smoothing algorithm is needed to obtain the smoothed distributions of ${\mathit{X}}_{0:L}$ .

Filtering algorithm

For a nonlinear system state estimation problem, combining GPSSMs and Bayesian filtering algorithms, a variety of state estimation algorithms are proposed. Ko and Fox ¹⁶ present a unified framework including the extended Kalman filter (EKF), the unscented Kalman filter (UKF), and the particle filter (PF). And, Deisenroth et al. ¹⁷ propose the Gauss process–assumed density filter (GP-ADF) to acquire the state estimation mean and covariance matrix accurately.

Combining GP-BayesFilters with the Rauch–Tung–Striebel (RTS) smoother, the corresponding RTS smoothing algorithms can be obtained directly. These algorithms use the GP based on the noise-free test input, and the predicted covariance matrix is a diagonal matrix only, which cannot describe the correlation among the output data. Deisenroth et al. ¹⁹ propose an analytical RTS smoother called as GP-RTSS algorithm, which is able to predict based on the noise test input and is proved to be more robust. The GP-RTSS algorithm is obtained by combining the RTSS algorithm with the GP-ADF algorithm.

The filtered distribution of the system state ${\mathit{x}}_k$ obeys the Gaussian distribution approximately, that is

$$p\left({\mathit{x}}_k|{\mathit{Y}}_{1:k}\right)\approx N({\mathit{x}}_k|{\hat{\mathit{x}}}_k,{\mathit{P}}_k)$$

(32)

$${\hat{\mathit{x}}}_k={\overline{\mathit{x}}}_k+{\mathit{P}}_k^{\mathrm{xy}}{\mathit{R}}_k^{-1}({\mathit{y}}_k-{\overline{\mathit{y}}}_k)$$

(33)

$${\mathit{P}}_k={\mathit{Q}}_k-{\mathit{P}}_k^{\mathrm{xy}}{\mathit{R}}_k^{-1}{\left({\mathit{P}}_k^{\mathrm{xy}}\right)}^T$$

(34)

where ${\hat{\mathit{x}}}_k$ is a filtered estimation mean vector, ${\mathit{P}}_k$ is a filtered estimation covariance matrix, ${\overline{\mathit{x}}}_k$ is a state prediction mean vector, ${\mathit{Q}}_k$ is a prediction state covariance matrix, ${\overline{\mathit{y}}}_k$ is a prediction measurement mean vector, ${\mathit{R}}_k$ is a prediction measurement covariance matrix, and ${\mathit{P}}_k^{\mathrm{xy}}$ is a cross-covariance matrix for the prediction state and the prediction measurement.

The filtered distribution of ${\mathit{x}}_{k-1}$ at the $(k-1)\mathrm{th}$ time is assumed as a priori Gaussian distribution, $N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})$ . ${\mathit{x}}_{k-1}$ as the noise is input into the state transition equation of GPSSM and then the prediction distribution of ${\mathit{x}}_k$ can be obtained as $N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)$ . If $N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)$ is a priori Gaussian distribution of ${\mathit{x}}_k$ , the prediction distribution of ${\mathit{y}}_k$ denoted as $N({\mathit{y}}_k|{\overline{\mathit{y}}}_k,{\mathit{R}}_k)$ can be obtained by inputting ${\mathit{x}}_k$ into the measurement equation of GPSSM. ${\mathit{P}}_k^{\mathrm{xy}}$ in formulas (33) and (34) is derived as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_k^{\mathrm{xy}}=\int \int ({\mathit{x}}_k-{\overline{\mathit{x}}}_k){({\mathit{y}}_k-{\overline{\mathit{y}}}_k)}^{T}p\left({\mathit{x}}_k,{\mathit{y}}_k|{\mathit{Y}}_{1:k-1}\right)d{\mathit{x}}_{k}d{\mathit{y}}_k}\\ \multicolumn{1}{c}{=\int \int {\mathit{x}}_k{\mathit{y}}_k^{T}p\left({\mathit{x}}_k,{\mathit{y}}_k|{\mathit{Y}}_{1:k-1}\right)d{\mathit{x}}_{k}d{\mathit{y}}_k-{\overline{\mathit{x}}}_k{\overline{\mathit{y}}}_k^T}\\ \multicolumn{1}{c}{=\int \int {\mathit{x}}_k{\mathit{y}}_k^{T}p\left({\mathit{y}}_k|{\mathit{x}}_k\right)p\left({\mathit{x}}_k|{\mathit{Y}}_{1:k-1}\right)d{\mathit{x}}_{k}d{\mathit{y}}_k-{\overline{\mathit{x}}}_k{\overline{\mathit{y}}}_k^T}\\ \multicolumn{1}{c}{{=}^{\left(4\right)}\int {\mathit{x}}_k[\int {\mathit{y}}_k^T\underset{i=1}{\overset{n}{\Pi }}N({\mathit{y}}_k^{(i)}|G{P}_{\rho }({\mathit{x}}_k,D_{h_i}^{*}),G{P}_{\Sigma }({\mathit{x}}_k,D_{h_i}^{*})+{\sigma }_{\mathrm{yi}}^2)d{\mathit{y}}_k]}\\ \multicolumn{1}{c}{N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k-{\overline{\mathit{x}}}_k{\overline{\mathit{y}}}_k^T}\\ \multicolumn{1}{c}{=\int {\mathit{x}}_k(G{P}_{\rho }({\mathit{x}}_k,D_{h_1}^{*}),G{P}_{\rho }({\mathit{x}}_k,D_{h_2}^{*}),\dots ,G{P}_{\rho }({\mathit{x}}_k,D_{h_m}^{*}))}\\ \multicolumn{1}{c}{N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k-{\overline{\mathit{x}}}_k{\overline{\mathit{y}}}_k^T}\\ \multicolumn{1}{c}{=[\int {\mathit{x}}_k{\mathrm{GP}}_{\rho }({\mathit{x}}_k,D_{h_1}^{*})N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k,\dots ,\int {\mathit{x}}_k{\mathrm{GP}}_{\rho }({\mathit{x}}_k,D_{h_m}^{*})}\\ \multicolumn{1}{c}{N({\mathit{x}}_k|{\overline{x}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k]-{\overline{x}}_k{\overline{\mathit{y}}}_k^T}\end{array}$$

(35)

In the step (4) of the above formula (35), $p({\mathit{y}}_k|{\mathit{x}}_k)$ is calculated by formula (21) and $p({\mathit{x}}_k|{\mathit{Y}}_{1:k-1})$ is the one-step prediction state distribution, that is $N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)$ . The key to calculate ${\mathit{P}}_k^{\mathrm{xy}}$ is the integral as follows, $i=1,2,\dots ,m$

$$\begin{array}{c}\multicolumn{1}{c}{\int {\mathit{x}}_{k}G{P}_{\rho }({\mathit{x}}_k,D_{h_i}^{*})N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{{=}^{\left(1\right)}\int {\mathit{x}}_k{\mathit{k}}_{h_i}^T{[{\mathit{K}}_{h_i}({\mathit{X}}_{1:L}^{*},{\mathit{X}}_{1:L}^{*})+{\sigma }_{\mathrm{yi}}^2\mathit{I}]}^{-1}}\\ \multicolumn{1}{c}{{\mathit{y}}_{1:L}^{(i)}N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{{=}^{\left(2\right)}\int {\mathit{x}}_k{\mathit{k}}_{h_i}^T\mathbf{\lambda }\left(D_{h_i}^{*}\right)N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{{=}^{\left(3\right)}\sum _{j=1}^L{\mathbf{\lambda }}^{(j)}\left(D_{h_i}^{*}\right)\int {\mathit{x}}_k{k}_{h_i}({\mathit{x}}_k,{\mathit{x}}_j^{*})N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{{=}^{(4)}\sum _{j=1}^L{c}_{\mathrm{ij}}{\mathbf{\lambda }}^{(j)}\left(D_{h_i}^{*}\right){({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*})}\end{array}$$

(36)

In step (1) of formula (36), $G{P}_{\rho }({\mathit{x}}_k,D_{h_i}^{*})$ can be obtained by replacing ${\mathit{\alpha }}_{*}$ and $〈\mathit{A},\mathit{\beta }〉$ with ${\mathit{x}}_k$ and $D_{h_i}^{*}=〈{\mathit{X}}_{1:L}^{*},{\mathit{y}}_{1:L}^{*(i)}〉$ in formula (4), respectively. In step (2) of formula (36), let $\mathbf{\lambda }(D_{h_i}^{\ast })={[{\mathit{K}}_{h_i}\left(X_{1:L}^{\ast },{\mathit{X}}_{1:L}^{\ast }\right)+{\sigma }_{yi}^2\mathit{I}]}^{-1}{\mathit{y}}_{1:L}^{\ast (i)}$ . In the step (3), let ${\mathit{k}}_{h_i}^T=(k_{h_i}({\mathit{x}}_k,{\mathit{x}}_1^{*}),k_{h_i}({\mathit{x}}_k,{\mathit{x}}_2^{*}),\dots ,k_{h_i}({\mathit{x}}_k,{\mathit{x}}_L^{*}))$ , where $k_{h_i}({\mathit{x}}_k,{\mathit{x}}_j^{*})$ is described by formula (16). Step (3) to step (4) in formula (36) is deducted as follows

$$\begin{array}{c}\multicolumn{1}{c}{\int {\mathit{x}}_k{k}_{h_i}({\mathit{x}}_k,{\mathit{x}}_j^{*})N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{=\int {\mathit{x}}_k{\sigma }_{h_i}^2\exp (-\frac{1}{2}{({\mathit{x}}_k-{\mathit{x}}_j^{*})}^T{\mathit{W}}_{h_i}^{-1}({\mathit{x}}_k-{\mathit{x}}_j^{*}))\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|{\mathit{Q}}_k\right|}^{\frac{1}{2}}}\exp (-\frac{1}{2}{({\mathit{x}}_k-{\overline{\mathit{x}}}_k)}^T{\mathit{Q}}_k^{-1}({\mathit{x}}_k-{\overline{\mathit{x}}}_k))d{\mathit{x}}_k}\\ \multicolumn{1}{c}{=\frac{{\sigma }_{h_i}^2}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|{\mathit{Q}}_k\right|}^{\frac{1}{2}}}\int {\mathit{x}}_k\exp (-\frac{1}{2}{({\mathit{x}}_k-{\mathit{x}}_j^{*})}^T{\mathit{W}}_{h_i}^{-1}({\mathit{x}}_k-{\mathit{x}}_j^{*})-\frac{1}{2}{({\mathit{x}}_k-{\overline{\mathit{x}}}_k)}^T{\mathit{Q}}_k^{-1}({\mathit{x}}_k-{\overline{\mathit{x}}}_k))d{\mathit{x}}_k}\\ \multicolumn{1}{c}{=\frac{{\sigma }_{h_i}^2}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|{\mathit{Q}}_k\right|}^{\frac{1}{2}}}\exp (-\frac{1}{2}{({\overline{\mathit{x}}}_k-{\mathit{x}}_j^{*})}^T{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\overline{\mathit{x}}}_k-{\mathit{x}}_j^{*}))}\\ \multicolumn{1}{c}{·\int {\mathit{x}}_k\exp \left(\begin{array}{c}-\frac{1}{2}{({\mathit{x}}_k-{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*}))}^T·{\left({\mathit{W}}_{h_i}{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}{\mathit{Q}}_k\right)}^{-1}\hfill \\ \multicolumn{1}{c}{·({\mathit{x}}_k-{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*}))}\end{array}\right)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{{=}^{\left(4\right)}{c}_{\mathrm{ij}}\int {\mathit{x}}_k\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|{\mathit{W}}_{h_i}{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}{\mathit{Q}}_k\right|}^{\frac{1}{2}}}}\\ \multicolumn{1}{c}{·\exp \left(\begin{array}{c}-\frac{1}{2}{({\mathit{x}}_k-{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*}))}^T·{\left({\mathit{W}}_{h_i}{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}{\mathit{Q}}_k\right)}^{-1}\hfill \\ \multicolumn{1}{c}{·({\mathit{x}}_k-{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*}))}\end{array}\right)d{\mathit{x}}_k}\\ \multicolumn{1}{c}{=c_{\mathrm{ij}}{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*})}\end{array}$$

(37)

$$\begin{array}{c}\multicolumn{1}{c}{c_{\mathrm{ij}}=\frac{{\sigma }_{h_i}^2{\left|{\mathit{W}}_{h_i}{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}{\mathit{Q}}_k\right|}^{1/2}}{{\left|{\mathit{Q}}_k\right|}^{1/2}}\exp }\\ \multicolumn{1}{c}{(-\frac{1}{2}{({\overline{\mathit{x}}}_k-{\mathit{x}}_j^{*})}^T{({\mathit{W}}_{h_i}+{\mathit{Q}}_k)}^{-1}({\overline{\mathit{x}}}_k-{\mathit{x}}_j^{*}))}\end{array}$$

(38)

In step (3) of formula (37), the integral is the expectation operation on ${\mathit{x}}_k$ essentially, and then, the result is the mean vector $({\mathit{W}}_{h_i}+{\mathit{Q}}_k{)}^{-1}({\mathit{W}}_{h_i}{\overline{\mathit{x}}}_k+{\mathit{Q}}_k{\mathit{x}}_j^{*})$ . $c_{\mathrm{ij}}$ is a normalization constant. Thus, the GP-ADF can be obtained as Algorithm 1.

Algorithm 1. GP-ADF

Step 1. Initialize $G({\theta }_{\mu },D^{*})$ and $N({\mathit{x}}_0|{\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ ; Step 2. Take $N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})$ as noise that is input to $G({\theta }_{\mu },D^{*})$ ; Step 3. Calculate $N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)$ ; Step 4. Take $N({\mathit{x}}_k|{\overline{\mathit{x}}}_k,{\mathit{Q}}_k)$ as noise that is input to $G({\theta }_{\mu },D^{*})$ ; Step 5. Calculate $N({\mathit{y}}_k|{\overline{\mathit{y}}}_k,{\mathit{R}}_k)$ ; Step 6. Calculate ${\mathit{P}}_k^{\mathrm{xy}}$ ; Step 7. Get ${\mathit{y}}_k$ ; Step 8. Calculate ${\hat{\mathit{x}}}_k={\overline{\mathit{x}}}_k+{\mathit{P}}_k^{\mathrm{xy}}{\mathit{R}}_k^{-1}({\mathit{y}}_k-{\overline{\mathit{y}}}_k)$ ; Step 9. Calculate ${\mathit{P}}_k={\mathit{Q}}_k-{\mathit{P}}_k^{\mathrm{xy}}{\mathit{R}}_k^{-1}{\left({\mathit{P}}_k^{\mathrm{xy}}\right)}^T$ ; Step 10. let $k=k+1$ and return to (2).

Smoothing algorithm

Algorithm 1 accurately calculates the mean vector and the covariance matrix of the filtered estimation. Combining the GP-ADF with the RTS smoother, the distribution of the smoothed estimation can be calculated as follows

$$p\left({\mathit{x}}_{k-1}|{\mathit{Y}}_{1:L}\right)\approx N({\mathit{x}}_{k-1}|{\tilde{\mathit{x}}}_{k-1},{\tilde{\mathit{P}}}_{k-1})$$

(39)

$${\tilde{\mathit{x}}}_{k-1}={\hat{\mathit{x}}}_{k-1}+{\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}({\tilde{\mathit{x}}}_k-{\overline{\mathit{x}}}_k)$$

(40)

$${\tilde{\mathit{P}}}_{k-1}={\mathit{P}}_k+{\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}{({\tilde{\mathit{P}}}_k-{\mathit{Q}}_k)}^T{\left({\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}\right)}^T$$

(41)

After executing Algorithm 1, the only unknown quantity is the matrix ${\mathit{P}}_{k-1,k}^{\mathrm{xx}}$ in the above formulas (39)–(41), which is calculated as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{k-1,k}^{\mathrm{xx}}=\int \int ({\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{k-1}){({\mathit{x}}_k-{\overline{\mathit{x}}}_k)}^{T}p\left({\mathit{x}}_{k-1},{\mathit{x}}_k|{\mathit{Y}}_{1:k-1}\right)d{\mathit{x}}_{k-1}d{\mathit{x}}_k}\\ \multicolumn{1}{c}{=\int \int {\mathit{x}}_{k-1}{\mathit{x}}_k^{T}p\left({\mathit{x}}_{k-1},{\mathit{x}}_k|{\mathit{Y}}_{1:k-1}\right)d{\mathit{x}}_{k-1}d{\mathit{x}}_k-{\hat{\mathit{x}}}_{k-1}{\overline{\mathit{x}}}_k^T}\\ \multicolumn{1}{c}{=\int {\mathit{x}}_{k-1}[\int {\mathit{x}}_k^T\underset{i=1}{\overset{n}{\Pi }}N({\mathit{x}}_k^{(i)}|G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_i}^{*}),G{P}_{\Sigma }({\mathit{x}}_{k-1},D_{f_i}^{*})+{\sigma }_{\mathrm{xi}}^2)d{\mathit{x}}_k]}\\ \multicolumn{1}{c}{·N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})d{\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{k-1}{\overline{\mathit{x}}}_k^T}\\ \multicolumn{1}{c}{=\int {\mathit{x}}_{k-1}(G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_1}^{*}),G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_2}^{*}),\dots ,G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_n}^{*}))}\\ \multicolumn{1}{c}{·N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})d{\mathit{x}}_{k-1}-{\hat{\mathit{x}}}_{k-1}{\overline{\mathit{x}}}_k^T}\\ \multicolumn{1}{c}{{=}^{\left(5\right)}[\dots ,\int {\mathit{x}}_{k-1}G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_i}^{*})N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})d{\mathit{x}}_{k-1},\dots ]-{\hat{\mathit{x}}}_{k-1}{\overline{\mathit{x}}}_k^T}\end{array}$$

(42)

In step (5) of formula (42), the calculation of the integral term is similar to formula (36)

$$\begin{array}{c}\int {\mathit{x}}_{k-1}G{P}_{\rho }({\mathit{x}}_{k-1},D_{f_i}^{*})N({\mathit{x}}_{k-1}|{\hat{\mathit{x}}}_{k-1},{\mathit{P}}_{k-1})d{\mathit{x}}_{k-1}\hfill \\ =\sum _{j=0}^{L-1}{d}_{\mathrm{ij}}{\mathbf{\chi }}^{(j)}\left(D_{f_i}^{*}\right){({\mathit{W}}_{f_i}+{\mathit{P}}_k)}^{-1}({\mathit{W}}_{f_i}{\hat{\mathit{x}}}_k+{\mathit{P}}_k{\mathit{x}}_j^{*})\hfill \end{array}$$

(43)

$$\mathbf{\chi }\left(D_{f_i}^{*}\right)={[{\mathit{K}}_{f_i}({\mathit{X}}_{0:L-1}^{*},{\mathit{X}}_{0:L-1}^{*})+{\sigma }_{\mathrm{xi}}^2\mathit{I}]}^{-1}{\mathit{x}}_{1:L}^{(i)}$$

(44)

$$\begin{array}{c}\multicolumn{1}{c}{d_{\mathrm{ij}}=\frac{{\sigma }_{f_i}^2{\left|{\mathit{W}}_{f_i}{({\mathit{W}}_{f_i}+{\mathit{P}}_k)}^{-1}{\mathit{P}}_k\right|}^{1/2}}{{\left|{\mathit{P}}_k\right|}^{1/2}}\exp }\\ \multicolumn{1}{c}{(-\frac{1}{2}{({\hat{\mathit{x}}}_k-{\mathit{x}}_j^{*})}^T{({\mathit{W}}_{f_i}+{\mathit{P}}_k)}^{-1}({\hat{\mathit{x}}}_k-{\mathit{x}}_j^{*}))}\end{array}$$

(45)

Thus, the GP-RTSS can be described as Algorithm 2.

Algorithm 2. GP-RTSS

Step 1. Initialize $G({\theta }_{\mu },D^{*})$ and $N({\mathit{x}}_0|{\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ ; Step 2. Perform the GP-ADF from $k=1$ to $k=L$ , and save the filtered distribution and the prediction distribution; Step 3. When $k=L$ , let ${\tilde{\mathit{x}}}_k={\hat{\mathit{x}}}_k$ and ${\tilde{\mathit{P}}}_k={\mathit{P}}_k$ ; Step 4. Calculate ${\mathit{P}}_{k-1,k}^{\mathrm{xx}}$ ; Step 5. Input the saved filtered distribution, prediction distribution and ${\mathit{P}}_{k-1,k}^{\mathrm{xx}}$ into (6) and (7); Step 6. Calculate ${\tilde{\mathit{x}}}_{k-1}={\hat{\mathit{x}}}_{k-1}+{\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}({\tilde{\mathit{x}}}_k-{\overline{\mathit{x}}}_k)$ ; Step 7. Calculate ${\tilde{\mathit{P}}}_{k-1}={\mathit{P}}_k+{\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}{({\tilde{\mathit{P}}}_k-{\mathit{Q}}_k)}^T{\left({\mathit{P}}_{k-1,k}^{\mathrm{xx}}{\mathit{Q}}_{k-1}^{-1}\right)}^T$ ; Step 8. Let $k=k-1$ and return to Step 5.

Maximization step

Algorithm 2 can obtain the smoothed distribution of the GPSSM. Next, the joint smoothed distribution $p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})$ is required to calculate the CEF

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathit{x}}_{k+1},{\mathit{x}}_k|{\mathit{Y}}_{1:L})=p({\mathit{x}}_k|{\mathit{x}}_{k+1},{\mathit{Y}}_{1:L})p({\mathit{x}}_{k+1}|{\mathit{Y}}_{1:L})}\\ \multicolumn{1}{c}{=p({\mathit{x}}_k|{\mathit{x}}_{k+1},{\mathit{Y}}_{1:k})p({\mathit{x}}_{k+1}|{\mathit{Y}}_{1:L})}\end{array}$$

(46)

where $p({\mathit{x}}_{k+1}|{\mathit{Y}}_{1:L})$ is the smoothed distribution of ${\mathit{x}}_{k+1}$ . The distribution $p({\mathit{x}}_k|{\mathit{x}}_{k+1},{\mathit{Y}}_{1:k})$ is calculated by $p({\mathit{x}}_k,{\mathit{x}}_{k+1}|{\mathit{Y}}_{1:k})$

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathit{x}}_k,{\mathit{x}}_{k+1}|{\mathit{Y}}_{1:k})=N}\\ \multicolumn{1}{c}{(\begin{array}{c}{\mathit{x}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{x}}_{k+1}}\end{array}|\left(\begin{array}{c}{\hat{\mathit{x}}}_k\hfill \\ \multicolumn{1}{c}{{\overline{\mathit{x}}}_{k+1}}\end{array}\right),\left(\begin{array}{cc}{\mathit{P}}_k\hfill & {\mathit{P}}_{k,k+1}^{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\left({\mathit{P}}_{k,k+1}^{\mathrm{xx}}\right)}^T}& {\mathit{P}}_{k+1}\hfill \end{array}\right))}\end{array}$$

(47)

Then

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathit{x}}_k|{\mathit{x}}_{k+1},{\mathit{Y}}_{1:k})=N}\\ \multicolumn{1}{c}{({\mathit{x}}_k|{\hat{\mathit{x}}}_k+{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}({\mathit{x}}_{k+1}-{\overline{\mathit{x}}}_{k+1}),{\mathit{P}}_k-{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}{({\mathit{P}}_{k,k+1}^{\mathrm{xx}})}^T)}\\ \multicolumn{1}{c}{=N({\mathit{x}}_k|({\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}){\mathit{x}}_{k+1}+({\hat{\mathit{x}}}_k-{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}{\overline{\mathit{x}}}_{k+1}),}\\ \multicolumn{1}{c}{{\mathit{P}}_k-{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}({\mathit{P}}_{k,k+1}^{\mathrm{xx}}{)}^T)}\\ \multicolumn{1}{c}{{=}^{(3)}N({\mathit{x}}_k|{\mathit{H}}_k{\mathit{x}}_{k+1}+{\mathit{U}}_k,{\mathbf{\Gamma }}_k)}\end{array}$$

(48)

In step (3) of formula (48), ${\mathit{H}}_k={\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}$ , ${\mathit{U}}_k={\hat{\mathit{x}}}_k-{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}{\overline{\mathit{x}}}_{k+1}$ , and ${\mathbf{\Gamma }}_k={\mathit{P}}_k-{\mathit{P}}_{k,k+1}^{\mathrm{xx}}{\mathit{Q}}_{k+1}^{-1}({\mathit{P}}_{k,k+1}^{\mathrm{xx}}{)}^T$ . By the property of Gauss distribution

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathit{x}}_{k+1},{\mathit{x}}_k|{\mathit{Y}}_{1:L})=N}\\ \multicolumn{1}{c}{(\begin{array}{c}{\mathit{x}}_{k+1}\hfill \\ \multicolumn{1}{c}{{\mathit{x}}_k}\end{array}|\left(\begin{array}{c}{\tilde{\mathit{x}}}_{k+1}\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_k{\tilde{\mathit{x}}}_{k+1}+{\mathit{U}}_k}\end{array}\right),\left(\begin{array}{cc}{\tilde{\mathit{P}}}_k\hfill & {\tilde{\mathit{P}}}_k{\mathit{H}}_k^T\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_k{\tilde{\mathit{P}}}_k}& {\mathit{H}}_k{\tilde{\mathit{P}}}_k{\mathit{H}}_k^T+{\mathbf{\Gamma }}_k\hfill \end{array}\right))}\\ \multicolumn{1}{c}{=N(\begin{array}{c}{\mathit{x}}_{k+1}\hfill \\ \multicolumn{1}{c}{{\mathit{x}}_k}\end{array}|\left(\begin{array}{c}{\tilde{\mathit{x}}}_{k+1}\hfill \\ \multicolumn{1}{c}{{\tilde{\mathit{x}}}_k}\end{array}\right),\left(\begin{array}{cc}{\tilde{\mathit{P}}}_k\hfill & {\tilde{\mathit{P}}}_k{\mathit{H}}_k^T\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_k{\tilde{\mathit{P}}}_k}& {\tilde{\mathit{P}}}_{k+1}\hfill \end{array}\right))}\end{array}$$

(49)

Thus, $l_0$ , $l_1$ , and $l_2$ are expressed as

$$l_0=\int p_{{\theta }_{\mu }}({\mathit{x}}_0|{\mathit{Y}}_{1:L})\log N({\mathit{x}}_0|{\hat{\mathit{x}}}_0,{\mathit{P}}_0)d{\mathit{x}}_0$$

(50)

$$\begin{array}{c}\multicolumn{1}{c}{l_1=\sum _{i=0}^{L-1}\int \int p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})\log \underset{j=1}{\overset{n}{\Pi }}N({\mathit{x}}_{i+1}^{(j)}|G{P}_{\rho }({\mathit{x}}_i,D_{f_j}),G{P}_{\Sigma }({\mathit{x}}_i,D_{f_j})+{\sigma }_{\mathrm{xj}}^2)d{\mathit{x}}_{i}d{\mathit{x}}_{i+1}}\\ \multicolumn{1}{c}{=-\frac{1}{2}\sum _{i=0}^{L-1}\sum _{j=1}^n\int \int p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})(\frac{{({\mathit{x}}_{i+1}^{(j)}-G{P}_{\rho }({\mathit{x}}_i,D_{f_j}))}^2}{G{P}_{\Sigma }({\mathit{x}}_i,D_{f_j})+{\sigma }_{\mathrm{xj}}^2}-\frac{1}{2}\log (G{P}_{\Sigma }({\mathit{x}}_i,D_{f_j})+{\sigma }_{\mathrm{xj}}^2))d{\mathit{x}}_{i}d{\mathit{x}}_{i+1}+\mathrm{const}}\end{array}$$

(51)

$$\begin{array}{c}\multicolumn{1}{c}{l_2=\sum _{i=1}^L\int p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})\log \underset{j=1}{\overset{m}{\Pi }}N({\mathit{y}}_i^{(j)}|G{P}_{\rho }({\mathit{x}}_i,D_{h_j}),G{P}_{\Sigma }({\mathit{x}}_i,D_{h_j})+{\sigma }_{\mathrm{yj}}^2)d{\mathit{x}}_i}\\ \multicolumn{1}{c}{=\sum _{i=1}^L\sum _{j=1}^m\int p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})\log N({\mathit{y}}_i^{(j)}|G{P}_{\rho }({\mathit{x}}_i,D_{h_j}),G{P}_{\Sigma }({\mathit{x}}_i,D_{h_j})+{\sigma }_{\mathrm{yj}}^2)d{\mathit{x}}_i}\\ \multicolumn{1}{c}{=-\frac{1}{2}\sum _{i=1}^L\sum _{j=1}^m\int p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})(\frac{{({\mathit{y}}_i^{(j)}-G{P}_{\rho }({\mathit{x}}_i,D_{h_j}))}^2}{G{P}_{\Sigma }({\mathit{x}}_i,D_{h_j})+{\sigma }_{\mathrm{yj}}^2}-\frac{1}{2}\log (G{P}_{\Sigma }({\mathit{x}}_i,D_{h_j})+{\sigma }_{\mathrm{yj}}^2))d{\mathit{x}}_i+\mathrm{const}}\end{array}$$

(52)

The variables $l_1$ and $l_2$ cannot be calculated analytically, because they contain the unprocessed integral in the closed form. The MC numerical integral method can provide any approximation accuracy for the integral operation. For $l_1$ , the MC method is used to sample N particles $\{{\mathit{x}}_{i+1}^g,{\mathit{x}}_i^g{\}}_{g=1}^N$ from $p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})$ randomly. Then, $l_1$ is approximated by ${\hat{l}}_1$

$${\hat{l}}_1=-\frac{1}{2N}\sum _{g=1}^N\sum _{i=0}^{L-1}\sum _{j=1}^n(\frac{{({\mathit{x}}_{i+1}^{g(j)}-G{P}_{\rho }({\mathit{x}}_i^g,D_{f_j}))}^2}{G{P}_{\Sigma }({\mathit{x}}_i^g,D_{f_j})+{\sigma }_{\mathrm{xj}}^2}+\log (G{P}_{\Sigma }({\mathit{x}}_i^g,D_{f_j})+{\sigma }_{\mathrm{xj}}^2))+\mathrm{const}$$

(53)

Similarly, for $l_2$ , $\{{\mathit{x}}_i^g{\}}_{g=1}^N$ are sampled from $p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})$ . $l_2$ is approximated by ${\hat{l}}_2$

$${\hat{l}}_2=-\frac{1}{2N}\sum _g^N\sum _{i=1}^L\sum _{j=1}^m(\frac{{({\mathit{y}}_i^{(j)}-G{P}_{\rho }({\mathit{x}}_i^g,D_{h_j}))}^2}{G{P}_{\Sigma }({\mathit{x}}_i^g,D_{h_j})+{\sigma }_{\mathrm{yj}}^2}+\log (G{P}_{\Sigma }({\mathit{x}}_i^g,D_{h_j})+{\sigma }_{\mathrm{yj}}^2))+\mathrm{const}$$

(54)

Based on the above content, $l(\theta ,{\theta }_{\mu })$ can be approximately expressed as $\hat{l}(\theta ,{\theta }_{\mu })=l_0+{\hat{l}}_1+{\hat{l}}_2$ .

After $\hat{l}(\theta ,{\theta }_{\mu })$ is obtained, the work shifts to maximize it about $\theta$ in order to obtain a new iteration ${\theta }_{\mu +1}$ . $\theta$ consists of three parts, namely, $({\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ , $\{{\sigma }_{f_i},{\mathit{W}}_{f_j},{\sigma }_{\mathrm{xi}}{\}}_{i=1}^n$ and $\{{\sigma }_{h_i},{\mathit{W}}_{h_j},{\sigma }_{\mathrm{yi}}{\}}_{i=1}^m$ . $({\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ only appear in $l_0$ , $\{{\sigma }_{f_i},{\mathit{W}}_{f_j},{\sigma }_{\mathrm{xi}}{\}}_{i=1}^n$ only appear in ${\hat{l}}_1$ , and $\{{\sigma }_{h_i},{\mathit{W}}_{h_j},{\sigma }_{\mathrm{yi}}{\}}_{i=1}^m$ only appear in ${\hat{l}}_2$ . $\{{\sigma }_{f_i},{\mathit{W}}_{f_j},{\sigma }_{\mathrm{xi}}{\}}_{i=1}^n$ and $\{{\sigma }_{h_i},{\mathit{W}}_{h_j},{\sigma }_{\mathrm{yi}}{\}}_{i=1}^m$ are recorded as ${\theta }_f$ and ${\theta }_h$ , respectively. Then, the maximization of $\hat{l}(\theta ,{\theta }_{\mu })$ can be decomposed into maximizing $l_0$ , ${\hat{l}}_1$ , and ${\hat{l}}_2$ on $\{({\hat{\mathit{x}}}_0,{\mathit{P}}_0),{\theta }_f,{\theta }_h\}$ , respectively.

From the lemma of Deisenroth et al., ¹⁹ when ${\hat{\mathit{x}}}_0={\tilde{\mathit{x}}}_0$ and ${\mathit{P}}_0={\tilde{\mathit{P}}}_0$ , $l_0$ takes the maximum value, where ${\tilde{\mathit{x}}}_0$ and ${\tilde{\mathit{P}}}_0$ are the mean and covariance of the initial state smoothed estimation distribution, respectively. The partial derivatives of ${\hat{l}}_1$ and ${\hat{l}}_2$ with respect to ${\theta }_f$ and ${\theta }_h$ , respectively can be computed analytically, which allows for gradient-based parameter optimization. In this article, the variable metric method ¹⁴ is used to search the maximum. The correction matrices constructed by the variable metric method are all symmetric positive definite matrices, which guarantee that the search direction are all descending direction and each iteration makes the function value decrease. In fact, ${\theta }_{\mu +1}$ in the next iteration is not necessarily globally optimal, which only needs to meet the inequality $\hat{l}({\theta }_{\mu +1},{\theta }_{\mu })>\hat{l}({\theta }_{\mu },{\theta }_{\mu })$ that ensures $p({\mathit{Y}}_{1:L}|G({\theta }_{\mu +1},D^{*}))$ $>p({\mathit{Y}}_{1:L}|G({\theta }_{\mu },D^{*}))$ .

Final optimization algorithm

The final optimization algorithm is described formally in this section. It is called as the EM-GP-RTSS algorithm, because it mainly consists of the EM algorithm and the GP-RTSS algorithm.

The key of Algorithm 3 is to replace the difficult-to-maximize measurement likelihood function with the CEF, as the MLE of hyperparameters for the expected GPSSM can be acquired asymptotically through iterative optimizing the CEF. Meanwhile, Algorithm 3 outputs the smoothed estimation at each EM iteration. Then, Algorithm 3 can also be called as a joint state estimation and model optimization algorithm for the GPSSM.

Algorithm 3. EM-GP-RTSS

Step 1. Initialize ${\theta }_{\mu }={\theta }^{*}$ , $\mu =0$ ; Step 2. E-step: (a) Run Algorithm 1 and Algorithm 2 to obtain $p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})$ , for $i=0,2,\dots ,L-1$ ; (b) Calculate t $p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})$ for $i=0,2,\dots ,L-1$ by formula (46); Step 3. M-step: (a) Sample ${\{{\mathit{x}}_{i+1}^g,{\mathit{x}}_i^g\}}_{g=1}^N$ randomly from $p_{{\theta }_{\mu }}({\mathit{x}}_{i+1},{\mathit{x}}_i|{\mathit{Y}}_{1:L})$ to calculate ${\hat{l}}_1$ by formula (53); (b) Sample ${\left\{{\mathit{x}}_i^g\right\}}_{g=1}^N$ randomly from $p_{{\theta }_{\mu }}({\mathit{x}}_i|{\mathit{Y}}_{1:L})$ to calculate ${\hat{l}}_2$ by formula (54); (c) Calculate the gradients of ${\hat{l}}_1$ and ${\hat{l}}_2$ on ${\theta }_f$ and ${\theta }_h$ , respectively; (d) Apply the variable metric method to search ${\theta }_f^{*}$ and ${\theta }_h^{*}$ ; (e) Let ${\theta }_{\mu +1}=\{({\hat{\mathit{x}}}_0,{\mathit{P}}_0),{\theta }_f^{*},{\theta }_h^{*}\}$ ; Step 4. Judge the convergence. If $\hat{l}({\theta }_{\mu +1},{\theta }_{\mu })-\hat{l}({\theta }_{\mu },{\theta }_{\mu })>\epsilon (\epsilon >0)$ , $\mu =\mu +1$ and return to step 2; otherwise, end.

The computational complexity of Algorithm 3 depends on the two parts of the suboperation: the first is to compute the smoothed distribution; the second is to compute the approximate CEF. For the GPSSM described in this article, the computational complexity of GP-RTSS algorithm is $O(L{n}^3+2L{m}^3)$ , where L is the sampling number. The computational complexity scales linearly with the sampling number L. The approximate CEF is calculated by formulas (53) and (54), and the corresponding computation complexity is $O(\mathrm{NL}(n+m))$ , where N is the number of particles in the MC numerical integral. The threshold value $\epsilon$ is not easy to choose in practice, and the alternative method must set the maximum number of EM iterations. Such a method can often be used to balance the amount of calculation and estimation accuracy.

Simulation experiment

In this section, the computer simulation experiments are carried out to verify the performance of the EM-GP-RTSS algorithm. The noise pendulum system is selected as the simulated object by the GPSSM in these experiments. Because both its state transition equation and measurement equation are typical nonlinear forms, it is often used to verify the performance of the nonlinear filtering and smoothing algorithms.^17,19

In this study, all experiments are performed using MATLAB R2013a on an Intel i5 quad-core 2.90-GHz 64-bit machine with 4 GB RAM.

Experimental settings

For a noise pendulum system with the unit length and the unit mass, its standard form of the discrete-time nonlinear SSM can be established as follows

$$\begin{array}{c}\multicolumn{1}{c}{\underset{{\mathit{x}}_k}{\underbrace{\left(\begin{array}{c}{x}_{1,k}\hfill \\ \multicolumn{1}{c}{x_{2,k}}\end{array}\right)}}=\underset{\mathit{f}\left({\mathit{x}}_k\right)}{\underbrace{\left(\begin{array}{c}{x}_{1,k-1}+x_{2,k-1}\mathrm{\Delta }t\hfill \\ \multicolumn{1}{c}{x_{2,k-1}-g\sin \left(x_{1,k-1}\right)\mathrm{\Delta }t}\end{array}\right)}}+{\mathit{q}}_{k-1}}\\ \multicolumn{1}{c}{y_k=\underset{h\left({\mathit{x}}_k\right)}{\underbrace{\sin \left(x_{1,k}\right)}}+r_k}\end{array}$$

(55)

where $x_{1,k}$ is the angle at the $k\mathrm{th}$ time, $x_{2,k}$ is the angular velocity at the $k\mathrm{th}$ time, g is the gravitational acceleration, $y_k$ is the measurement at the $k-\mathrm{th}$ time, $\mathrm{\Delta }t$ is the sampling period of the sensor, ${\mathit{q}}_{k-1}~N({\mathit{q}}_{k-1}|0,\mathit{Q})$ , $r_k~N(r_k|0,R^2)$ and

$$\mathit{Q}=\left[\begin{array}{cc}\frac{q\mathrm{\Delta }{t}^3}{3}\hfill & \frac{q\mathrm{\Delta }{t}^2}{2}\hfill \\ \multicolumn{1}{c}{\frac{q\mathrm{\Delta }{t}^2}{2}}& q\mathrm{\Delta }t\hfill \end{array}\right]$$

(56)

where q is the spectral density.

On one hand, the nonlinear SSM (equation (55)) is served as a benchmark model to test the EM-GP-RTSS; on the other hand, it provides true states to train the GPSSM* that is taken as the GPSSM in the laboratory. We get different instantiated pendulum models (equation (55)) by setting different values for parameters $(q,R,g)$ . Five parameter combinations $(q,R,g)$ are set in the simulation experiment, as shown in Table 1.

Table 1.

Parameter settings of five pendulum models.

Parameters	q	R	g
Pendulum model0	0.01	0.01	9.831
Pendulum model1	0.05	0.05	9.831
Pendulum model2	0.1	0.1	9.831
Pendulum model3	0.5	0.5	9.831
Pendulum model4	0.01	0.01	9.719

Five pendulum models in Table 1 represent five homogeneous nonlinear dynamic systems. The pendulum model⁰ is used as the model in the laboratory that provides true states. The remaining four pendulum models describe noise pendulum systems in the different non-laboratory environments. For the pendulum models^1–3, we set the different values of $(q,R)$ from that of the pendulum model⁰, indicating the different noise from the laboratory. For the pendulum model⁴, a different gravitational acceleration g is set to represent a different geographical location from the laboratory, because g is a function of the latitude and the elevation. For all of the above pendulum models, set the sensor sampling period $\mathrm{\Delta }t=0.1\mathrm{s}$ , the initial time $t_0=0\mathrm{s}$ , and the sample number $L=50$ . The initial state obeys ${\mathit{x}}_0~N({\mathit{x}}_0|{\hat{\mathit{x}}}_0,{\mathit{P}}_0)$ , where ${\hat{\mathit{x}}}_0~N({\hat{\mathit{x}}}_0|0,{\mathit{E}}_{2\times 2})$ , ${\mathit{P}}_0=0.1{\mathit{E}}_{2\times 2}$ , ${\mathit{E}}_{2\times 2}=[1,0;0,1]$ .

By running the pendulum model⁰, $({\mathit{X}}_{0:50}^{*},{\mathit{Y}}_{1:50}^{*})$ are obtained, which are true system states and measurements. Then, the GPSSM* is trained with $({\mathit{X}}_{0:50}^{*},{\mathit{Y}}_{1:50}^{*})$ to simulate the pendulum model⁰. Trained hyperparameters of the GPSSM* are $\{{\sigma }_{f_i},{\mathit{W}}_{f_i},{\sigma }_{\mathrm{xi}}{\}}_{i=1}^2$ and $\{{\sigma }_h,{\mathit{W}}_h,{\sigma }_y\}$ where ${\mathit{W}}_{f_i}=[s_{f_{i,1}}^2,0;0,s_{f_{i,2}}^2]$ and ${\mathit{W}}_h=[s_{h_1}^2,0;0,s_{h_2}^2]$ . $\{{\mathit{X}}_{0:50}^i,{\mathit{Y}}_{1:50}^i{\}}_{i=1}^4$ are obtained by running the pendulum models^1–4 separately, where ${\mathit{X}}_{0:50}^i$ is unknown in the non-laboratory environment. If the GPSSM* is used to estimate ${\mathit{X}}_{0:50}^i$ given ${\mathit{Y}}_{1:50}^i$ , there is an overfitting issue. If a new GPSSMⁱ, which is obtained by optimizing the GPSSM* with the EM-GP-RTSS algorithm for the pendulum modelⁱ, can improve the estimation accuracy of ${\mathit{X}}_{0:50}^i$ , then it demonstrates that the GPSSMⁱ can simulate the pendulum modelⁱ well. For the pendulum models^1–4, this article uses the GPSSM*, the GPSSMⁱ, and the parametric model (equation (55)) to estimate their states. For the GPSSMs and the traditional SSMs, the Unscented Rauch–Tung–Striebel smoother (URTSS) ²⁰ that is a classical nonlinear smoothing algorithm is used as a benchmark for comparison with the GP-RTSS proposed in this article. The root mean square error (RMSE) is used to measure the accuracy of state estimation. The RMSE of the state $x_{1,k}$ is calculated as follows

$$\mathrm{RMSE}(x_{1,k})=\sqrt{\frac{1}{N}\sum _{i=1}^N{({\hat{x}}_{1,k}^i-x_{1,k})}^2}$$

(57)

where ${\hat{x}}_{1,k}^i$ is the smoothed estimation mean of $x_{1,k}$ at the $i\mathrm{th}$ MC simulation experiment. The RMSE of the state $x_{2,k}$ is similar to formula (57). In the EM-GP-RTSS, the convergence condition is set to 20 maximum iterations, and the number of sampling particles in the MC numerical integral is set as $N=200$ . The parameter setting of the URTSS algorithm is the same as Särkkä. ²⁰

Experimental results

True states ${\mathit{X}}_{0:50}^{*}$ of the pendulum model⁰ are shown in Figure 2. Measurements $\{{\mathit{Y}}_{1:50}^i{\}}_{i=0}^5$ of the pendulum models^0–4 are shown in Figure 3. $({\mathit{X}}_{0:50}^{*},{\mathit{Y}}_{1:50}^{*})$ are the training data of the GPSSM*. In simulation experiments, the BFGS algorithm ²¹ is used to train hyperparameters of the GPSSM*, which are recorded in Table 2. $({\mathit{X}}_{0:50}^{*},{\mathit{Y}}_{1:50}^{*})$ are the input data and the output data to construct the kernel functions of the GPSSMs^1–4. $\{{\mathit{Y}}_{1:50}^i{\}}_{i=1}^4$ are taken as measurements in these non-laboratory environments. The EM-GP-RTSS is adopted to optimize the GPSSM* in order to obtain the GPSSMs^1–4 whose hyperparameters are recorded in Table 3. The kernel functions $k_{f_1}$ , $k_{f_2}$ , and $k_h$ in Table 3 are defined in formulas (15) and (16).

Figure 2.

The true states obtained by the pendulum model⁰.

Figure 3.

The measurements obtained by the pendulum models^0–4.

Table 2.

Hyperparameters of GPSSM*.

Function	kf 1	kf 2	kh
${\sigma }_{f(h)}$	10.435	3.5634	1.7309
$s_{f{(h)}_1}$	52.625	14.844	95.455
$s_{f{(h)}_2}$	3.4048	7.5330	0.6710
${\sigma }_{x(y)}$	0.0035	0.0327	0.0086

GPSSM: Gauss process state-space model.

Table 3.

Hyperparameters of GPSSM^1–4.

Model	GPSSM1			GPSSM2			GPSSM3			GPSSM4
Function	kf 1	kf 2	kh	kf 1	kf 2	kh	kf 1	kf 2	kh	kf 1	kf 2	kh
${\sigma }_{f(h)}$	18.99	4.099	2.163	24.05	4.801	1.640	24.69	3.151	1.444	7.655	4.718	1.838
$s_{f{(h)}_1}$	104.7	14.42	66.80	96.64	34.39	64.68	106.7	5.535	39.65	71.29	40.76	100.2
$s_{f{(h)}_2}$	8.003	6.933	1.013	11.10	17.13	0.764	11.61	3.769	0.814	3.680	19.85	0.763
${\sigma }_{x(y)}$	0.010	0.074	0.052	0.015	0.106	0.098	0.025	0.243	0.466	0.004	0.037	0.009

GPSSM: Gauss process state-space model.

Figures 4–7 show the smoothed estimation mean of pendulum models^1–4 obtained after running the EM-GP-RTSS to optimize the GPSSM*, respectively. And, the estimation of the GP-RTSS with the GPSSM* and the estimation of the URTSS with the true SSM (53) are also given. From Figures 4–7, it can be seen that the smoothed estimation mean of the EM-GP-RTSS can more accurately follow true states for the pendulum models^1–4 than that of the GP-RTSS with the GPSSM*, which can also be seen from the corresponding standard deviation in Table 4. Therefore, it demonstrates that directly applying the GPSSM* to simulate these non-laboratory pendulum models^1–4 introduces the model bias. It can be seen from Table 4 that the GP-RTSS with GPSSMs^1–4 can reach or exceed the accuracy level of the URTSS with the SSMs^1–4.

Figure 4.

The smoothed estimation results for the pendulum model¹: (a) angle smoothing estimation mean and (b) angular velocity smoothing estimate mean.

Figure 5.

The smoothed estimation results for the pendulum model²: (a) angle smoothing estimation mean and (b) angular velocity smoothing estimate mean.

Figure 6.

The smoothed estimation results for the pendulum model³: (a) angle smoothing estimation mean and (b) angular velocity smoothing estimate mean.

Figure 7.

The smoothed estimation results for the pendulum model⁴: (a) angle smoothing estimation mean and (b) angular velocity smoothing estimate mean.

Table 4.

Standard deviation of smoothed estimation mean.

Algorithm	EM-GP-RTSS		GP-RTSS		URTSS
Standard deviation	Angle	Angular velocity	Angle	Angular velocity	Angle	Angular velocity
Pendulum model 1	0.0666	0.1340	0.0890	0.1381	0.1101	0.2352
Pendulum model 2	0.3547	0.3450	0.3996	0.4196	0.1434	0.2879
Pendulum model 3	0.4197	0.4299	0.6401	1.0124	0.8329	0.8476
Pendulum model 4	0.0239	0.0581	0.0447	0.1385	0.0430	0.0433

EM: expectation–maximization; GP: Gauss process; RTSS: Rauch–Tung–Striebel smoother; URTSS: unscented Rauch–Tung–Striebel smoother.

The one-step prediction performances of GPSSMs^1–4 and the GPSSM* are tested in order to further discuss the advantage of GPSSMs^1–4. Run pendulum models^1–4 again in the time period $[0,1\mathrm{s}]$ and generate four sets of new true states and measurements denoted as $\{{{\mathit{X}}^{\prime }}_{0:10}^i,{{\mathit{Y}}^{\prime }}_{1:10}^i{\}}_{i=1}^4$ . $\{{{\mathit{X}}^{\prime }}_{0:9}^i{\}}_{i=1}^4$ are input to the state transition equations of GPSSMs^1–4 and then the prediction distributions of $\{{{\mathit{X}}^{\prime }}_{1:10}^i{\}}_{i=1}^4$ are obtained from formula (20). Similarly, $\{{{\mathit{X}}^{\prime }}_{1:10}^i{\}}_{i=1}^4$ are input to the measurement equations of GPSSMs^1–4 and then the prediction distribution of $\{{{\mathit{Y}}^{\prime }}_{1:10}^i{\}}_{i=0}^4$ are obtained from formula (21).

The prediction distributions from the GPSSMs^1–4 are compared with the prediction distributions from the GPSSM*, as shown in Figures 8–11. The error bars show the twice standard deviations of the state prediction distributions or the measurement prediction distributions. The red circle indicates the corresponding true value of the prediction. The distributions obtained by the GPSSMs^1–4 can contain the true values within the twice standard deviations. The true value appears in the range of the twice standard deviations of the distributions from the GPSSMs^1–4 with the probability 95.5%. Except the pendulum model⁴, the prediction distributions of the GPSSM* cannot contain the true value in the twice standard deviations. When the noise of a nonlinear dynamic system increases in these non-laboratory environments, the confidence of the prediction distribution with the GPSSM* is lower than that with the GPSSMs^1–4 optimized by the EM-GP-RTSS. And, the mean values of the predicted distributions with the GPSSMs^1–4 are closer to the true value.

Figure 8.

The one-step prediction distributions with GPSSM* and GPSSM¹ for the pendulum model¹. (a) The prediction distribution of the state transition equation GP model. (b) The prediction distribution of the measurement equation GP model.

Figure 9.

The one-step prediction distributions with GPSSM* and GPSSM² for the pendulum model². (a) The prediction distribution of the state transition equation GP model. (b) The prediction distribution of the measurement equation GP model.

Figure 10.

The one-step prediction distributions with GPSSM* and GPSSM³ for the pendulum model³. (a) The prediction distribution of the state transition equation GP model. (b) The prediction distribution of the measurement equation GP model.

Figure 11.

The one-step prediction distribution with GPSSM* and GPSSM⁴ for the pendulum model⁴. (a) The prediction distribution of the state transition equation GP model. (b) The prediction distribution of the measurement equation GP model.

The GPSSMs^1–4 is reused for the state estimation problem of the corresponding nonlinear dynamic system to evaluate the improvement of simulation performance. Execute 100 MC simulations for the each GPSSM of GPSSMs^1–4. The true states and the measurements are generated by the pendulum models^1–4, respectively. The RMSE of angle and angular velocity are calculated. The URTSS with the traditional SSM provides a benchmark for comparison. For the GPSSMs, the GP-RTSS is used to calculate the smoothed estimation.

Figures 12–15 show the RMSE for the state estimation of the pendulum models^1–4 using the GPSSM*, the GPSSMs^1–4 and the traditional SSMs^1–4, respectively. It can be seen that the optimized GPSSMs^1–4 have the lower RMSEs than the GPSSM* in both state dimensions. It shows that the optimized GPSSMs^1–4 can better simulate the corresponding pendulum models^1–4 and reduce the model bias introduced by the GPSSM*. Meanwhile, the GP-RTSS with the GPSSMs can achieve or exceed the state estimation accuracy of the URTSS with the traditional SSMs.

Figure 12.

Comparison of state estimation accuracy with GPSSM*, GPSSM¹, and traditional SSM¹ for pendulum model¹ in 100 MC simulations: (a) RMSE of angle and (b) RMSE of angular velocity.

Figure 13.

Comparison of state estimation accuracy with GPSSM*, GPSSM², and traditional SSM² for pendulum model² in 100 MC simulations: (a) RMSE of angle and (b) RMSE of angular velocity.

Figure 14.

Comparison of state estimation accuracy with GPSSM*, GPSSM³, and traditional SSM³ for pendulum model³ in 100 MC simulations: (a) RMSE of angle and (b) RMSE of angular velocity.

Figure 15.

Comparison of state estimation accuracy with GPSSM*, GPSSM⁴, and traditional SSM⁴ for pendulum model⁴ in 100 MC simulations: (a) RMSE of angle and (b) RMSE of angular velocity.

Figure 16 shows the standard deviation of the smoothed estimation at each iteration of optimizing the GPSSM* by the EM-GP-RTSS to simulate the pendulum model³. When optimizing the GPSSM³ by the EM-GP-RTSS iteratively, the smoothed estimation can also be obtained, whose error constantly decreases and eventually converges.

Figure 16.

Standard deviation of smoothed estimation by the EM-GP-RTSS at each iteration: (a) standard deviation of angle and (b) standard deviation of angular velocity.

In the simulation experiment, the convergence condition of the EM-GP-RTSS is set to judge whether the maximum number of iterations is reached. Setting a loose maximum number of iterations can avoid the EM falling into a local optimum and directly control the amount of computation. Because of the operability, it is widely used to solve the MLE problem by the EM.^11,13

Conclusion and future works

In this article, for the nonlinear dynamic systems, a novel optimization algorithm for the GPSSM is proposed by combining the EM algorithm with the GP-RTSS algorithm, that is, the EM-GP-RTSS algorithm.

The main works are as follows: first, according to the concept of the GP, a theoretical formulation of the GPSSM is proposed, which is not found in previous references. The description provides a rigorous mathematical basis for the development of the EM-GP-RTSS algorithm. Second, according to the maximum likelihood estimation (MLE) criterion, a GPSSM optimization framework with the EM algorithm is proposed. In the EM algorithm, the unknown system state is considered as the lost data, and the maximization of measurement likelihood function is transformed into that of the CEF. Then, the GP-ADF algorithm and the GP-RTSS algorithm are proposed with the GPSSM defined in this article, in order to calculate the smoothed distribution in the CEF. Finally, the MC numerical integral method is used to obtain the approximate expression of the CEF.

In the simulation experiment, a noise pendulum model is used as a representative of the nonlinear dynamic system to verify the effectiveness of the EM-GP-RTSS algorithm proposed in this article. The noise pendulum model⁰ is selected as the nonlinear system in the laboratory, whose true states are obtained, and then a GPSSM* is trained to simulate it. The EM-GP-RTSS algorithm is used to optimize the GPSSM* for simulating four pendulum models^1–4 to get the GPSSMs^1–4. The experimental results demonstrate that the GPSSMs^1–4 can estimate the system state of the pendulum models^1–4 more accurately than the GPSSM* and can reach the estimation accuracy of the traditional SSM. The one-step prediction distributions from the GPSSMs^1–4 have the higher confidence than those of the GPSSM*. The EM-GP-RTSS algorithm is helpful to expand the application scope of the GPSSM from a laboratory environment to a practical non-laboratory environment.

In future research, first, a linear kernel function could be designed to develop a corresponding smoothing algorithm for the GPSSMs. A smoothing algorithm based on importance sampling methods can be used to obtain the nonlinear non-Gaussian smoothed estimation. Second, it is necessary to balance the computational complexity by selecting an appropriate number of particles in the MC numerical integral method. However, the relationship between the number of particles and the accuracy is difficult to determine. In addition to the MC numerical integration method, the fixed-point sampling numerical integration methods, such as based on the unscented transition or the spherical-radial cubature rule, may be adopted. The advantage of fixed-point sampling numerical integration methods is to require the fewer sampling points compared with the MC numerical integration method. Then, we will study how to use the high-degree fixed-point sampling numerical integration method to obtain the high-accuracy approximations for formulas (51) and (52).