# Data-Driven Identification of Noise Covariances in Kalman Filtering for Virtual Sensing Applications

Szymon Gres, Michael Döhler, Vasilis Dertimanis, Eleni Chatzi

► To cite this version:

Szymon Gres, Michael Döhler, Vasilis Dertimanis, Eleni Chatzi. Data-Driven Identification of Noise Covariances in Kalman Filtering for Virtual Sensing Applications. IOMAC 2024 - 10th International Operational Modal Analysis Conference, May 2024, Naples, Italy. pp.375-382, 10.1007/978-3-031-61425-5_36. hal-04639377

HAL Id: hal-04639377

https://inria.hal.science/hal-04639377v1

Submitted on 9 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Data-driven identification of noise covariances in Kalman filtering for virtual sensing applications

Szymon Gres<fnref n="1" />,<fnref n="3" />, Michael Döhler<fnref n="2" />, Vasilis Dertimanis<fnref n="3" />, and Eleni Chatzi<fnref n="3" />

**Abstract.** The optimality of the Kalman filter for state-estimation depends on the knowledge of the process and measurement noise covariance. In applications, these covariances are often treated as tuning parameters, often adjusted in a heuristic manner based on user-defined performance criteria. While several methods to identify them from data exist, some require the use of optimization algorithms, or inversion of large matrices, which is numerically inefficient. In this paper we review a simple data-driven subspace identification approach to estimate the process and measurement noise covariance, and apply it in the context of virtual sensing. The estimates of the noise covariances are obtained from the residuals of a regression of the banks of shifted Kalman filter states obtained from data and the model matrices of the monitored system. The performance of the approach is illustrated on a numerical example of a chain system by comparing the system states in unmeasured locations to the Kalman filter-based predictions.

## 1 Introduction

In state estimation the dynamic state of the system is obtained from process data and a system model. The process data constitutes input and output measurements at locations observed through sensors, which are a low-resolution subset of the full set of dynamic inputs and responses. The unmeasured inputs and the inherent errors associated with measurements give birth to the process and measurement noise that perturb the state and output equations, affecting state estimation. The Kalman filter is a well-suited tool for state estimation under the premise that the covariances of process and measurement noise are known. In applications, however, the noise covariance are generally unknown and their values are oftentimes tuned to optimize user-defined criteria that characterize filtering performance. While filter tuning is usually associated with heuristic optimization approaches, at its core lies a noise covariance identification problem.

Identification of the noise covariances has been addressed in numerous papers, and typically is divided into four categories of methods: Bayesian [5, 14], maximum likelihood [3, 15], covariance matching [8, 11] and correlation-based

<footnote n="1">¹ Department of Civil and Architectural Engineering, Aarhus University, 8000 Aarhus, Denmark</footnote>

<footnote n="2">² Univ. Gustave Eiffel, Inria, COSYS/SII, I4S, 35042 Rennes, France</footnote>

<footnote n="3">³ Department of Civil, Environmental and Geomatic Engineering, ETH Zürich, 8093 Zürich, Switzerland</footnote>
approaches. The correlation-based methods yield efficient estimates of the noise covariances and thus has found great interest in the engineering community. Autocovariance least-squares (ALS) [9, 1] is a promising correlation-based technique whose core is to estimate the covariance of the process and measurement noise from a covariance matrix of an innovation sequence. The estimation is formulated as a linear least-squares problem with additional constraints to ensure positive semidefiniteness of the estimated covariance matrices [9]. This optimization problem is convex and can be solved using e.g., Newton's method [9], or interior-point algorithms [1].

To circumvent using costly numerical optimization algorithms, the noise covariances can also be estimated from the residuals of state sequences identified with subspace-based identification methods [10]. The main drawback of the subspace-based approach is that the covariances are identified up to a similarity transform, thus cannot be directly used for a model-based Kalman filter for state estimation. To address this, the estimation basis is changed to the model basis with a linear regression in [7]. The resultant process and measurement covariance matrices are positive semi-definite and can readily be used for state estimation.

The purpose of this paper is to compare the subspace-based approach with the ALS method in a virtual sensing application. The paper is organized as follows. The system modelling is defined in Section 2, the subspace-based approach to estimate noise covariance is recalled in Section 3 and the numerical application is presented in Section 4.

## 2 Modelling

Consider the state-space system

$$x_{k+1} = Ax_k + Bu_k + w_k \quad (1)$$

$$y_k = Cx_k + Du_k + v_k, \quad (2)$$

where $x_k \in \mathbb{R}^n$ are the states, $y_k \in \mathbb{R}^r$ are the outputs, $u_k \in \mathbb{R}^u$ are the inputs, and $A \in \mathbb{R}^{n \times n}$, $C \in \mathbb{R}^{r \times n}$, $B \in \mathbb{R}^{n \times u}$, $D \in \mathbb{R}^{r \times u}$ are the state transition, observation, input and feedthrough matrices, respectively. The process noise $w_k \in \mathbb{R}^n$ and the measurement noise $v_k \in \mathbb{R}^r$ are modelled as zero-mean, Gaussian and white processes, of covariance and cross-covariance matrices

$$\begin{bmatrix} Q & S \\ S^T & R \end{bmatrix} = \mathbf{E} \left( \begin{bmatrix} w_k \\ v_k \end{bmatrix} \begin{bmatrix} w_k^T & v_k^T \end{bmatrix} \right). \quad (3)$$

The recursive Kalman predictor to obtain the one step-ahead state estimates $\hat{x}_{k+1|k}$ writes

$$\hat{x}_{k+1|k} = A\hat{x}_{k|k-1} + Bu_k + L_k (y_k - C\hat{x}_{k|k-1} - Du_k) \quad (4)$$

$$P_{k+1|k} = AP_{k|k-1}A^T + Q - L_k (AP_{k|k-1}C^T + S)^T, \quad (5)$$

$$L_k = (AP_{k|k-1}C^T + S)(CP_{k|k-1}C^T + R)^{-1}, \quad (6)$$
where $L_k$ is the predictor gain, and $P_{k+1|k} = E((\hat{x}_{k+1|k} - x_{k+1})(\hat{x}_{k+1|k} - x_{k+1})^T)$ is the prediction error covariance.

The goal of this paper is to obtain estimates of the triplet of the noise covariance matrices $(Q, R, S)$ from data-model residuals obtained through subspace identification.

### 3 Subspace-based noise covariance estimation

Usually, subspace-based system identification methods are used for the estimation of the system matrices $(A, B, C, D)$ from input/output measurements [12, 13]. Here, their properties will be used to obtain estimates $\hat{Q}, \hat{R}, \hat{S}$ from the measurements, where $(A, B, C, D)$ are already known from a model. This has been proposed in [7] based on the input/output N4SID subspace method [10].

**Definition 1 (Data matrix).** Let $a_k \in \mathbb{R}^b$ be a discrete signal at time step $k$. The parameter $p$ defines the 'past' and 'future' data horizons. For $0 \le i \le j \le 2p-1$ the data matrix $A_{i|j}$ writes

$$A_{i|j} = \frac{1}{\sqrt{N}} \begin{bmatrix} a_i & a_{i+1} & \dots & a_{i+N-1} \\ a_{i+1} & a_{i+2} & \dots & a_{i+N} \\ \vdots & \vdots & \ddots & \vdots \\ a_j & a_{j+1} & \dots & a_{j+N-1} \end{bmatrix} \in \mathbb{R}^{(j-i+1)b \times N}. \quad (7)$$

From the input and output data, define the data matrices $\mathcal{Y}^- \in \mathbb{R}^{pr \times N}$, $\mathcal{Y}^+ \in \mathbb{R}^{pr \times N}$, $\mathcal{U}^- \in \mathbb{R}^{pu \times N}$, $\mathcal{U}^+ \in \mathbb{R}^{pu \times N}$ and $\mathcal{W}^- \in \mathbb{R}^{p(u+r) \times N}$

$$\mathcal{Y}^{-} = \mathcal{Y}_{0|p-1}, \quad \mathcal{Y}^{+} = \mathcal{Y}_{p|2p-1}, \quad \mathcal{U}^{-} = \mathcal{U}_{0|p-1}, \quad \mathcal{U}^{+} = \mathcal{U}_{p|2p-1}, \quad \mathcal{W}^{-} = [\mathcal{U}^{-T} \mathcal{Y}^{-T}]^{T}. \quad (8)$$

Furthermore, denote respectively the future and past block-row matrices for the Kalman filter states as $\hat{\mathcal{X}}^{-} = \hat{\mathcal{X}}_{0|0}$, $\hat{\mathcal{X}}^{+} = \hat{\mathcal{X}}_{p|p}$. The observability matrix $\Gamma \in \mathbb{R}^{pr \times n}$

$$\Gamma = \begin{bmatrix} C \\ CA \\ CA^2 \\ \vdots \\ CA^{p-1} \end{bmatrix} \quad (9)$$

and banks of non steady-state Kalman filter state sequence $\hat{\mathcal{X}}$ can be obtained from a subspace of a projection of data matrices $\mathcal{Y}^{-}, \mathcal{Y}^{+}, \mathcal{U}^{-}, \mathcal{U}^{+}$ after [10].

The N4SID algorithm [10] is based on the projection of the input and output measurement matrices

$$\mathcal{H} = \mathcal{Y}^{+}/\mathcal{U}^{+}\mathcal{W}^{-} = (\mathcal{Y}^{+}/\mathcal{U}^{+\perp})(\mathcal{W}^{-}/\mathcal{U}^{+\perp})^{\dagger}\mathcal{W}^{-} \quad (10)$$
where the data matrices are defined in (8), (·)† denotes the pseudoinverse, and the orthogonal projection of a data matrix A onto the orthogonal complement of B is defined by A/B⊥ = A − A/B = A − AB$^{T}$(BB$^{T}$)$^{−1}$B. Matrix H factors into observability matrix and the Kalman state sequence

$$
\mathcal{H} \approx \Gamma \hat{\chi}^{+} \qquad (11)
$$

where the estimates can be retrieved from SVD

$$
\mathcal{H} = [U_1 \ U_2] \begin{bmatrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{bmatrix} \begin{bmatrix} V_1^T \\ V_2^T \end{bmatrix}, \quad \hat{\Gamma} = U_1 \Sigma_1^{1/2}, \quad \hat{\chi}^+ = \Sigma_1^{1/2} V_1^T. \quad (12)
$$

Furthermore, the Kalman state sequence shifted by one sample $\hat{\chi}_s^+ = \hat{\chi}_{p+1|p+1}$ can be obtained based on the shifted input and output measurement matrices $\mathcal{Y}_s^- = \mathcal{Y}_{0|p-2}$, $\mathcal{Y}_s^+ = \mathcal{Y}_{p-1|2p-1}$, $\mathcal{U}_s^- = \mathcal{U}_{0|p-2}$, $\mathcal{U}_s^+ = \mathcal{U}_{p-1|2p-1}$, $\mathcal{W}_s^- = [\mathcal{U}_s^{-T} \mathcal{Y}_s^{-T}]^T$ from the projection of the shifted matrices $\mathcal{H}_s = \mathcal{Y}_s^+ / \mathcal{U}_s^+ \mathcal{W}_s^-$ by

$$
\hat{\chi}_s^+ = \hat{\Gamma}_s^\dagger \mathcal{H}_s, \qquad (13)
$$

where $\hat{\Gamma}_s$ is the observability matrix $\hat{\Gamma}$ without the last block row.

The estimated Kalman state sequence $\hat{\chi}^+$ and its shifted version $\hat{\chi}_s^+$ are now used to infer about the noise covariances (3) corresponding to system (1)–(2) that is defined with the model-based system matrices $(A, B, C, D)$. However, in decomposition (12) the estimated state sequences correspond in general to a different state-space basis than the one defined by $(A, B, C, D)$. This can be easily fixed by estimating the change of basis matrix $\hat{T}$ that relates the estimated observability matrix $\hat{\Gamma}$ to the observability matrix $\Gamma$ defined from the model matrices $A$ and $C$ (see (9)) by $\Gamma\hat{T} = \hat{\Gamma}$, namely

$$
\hat{T} = \Gamma^{\dagger} \hat{\Gamma}. \qquad (14)
$$

Then, the residuals for the noise covariance estimation can be obtained based
on the shifted state sequences and the model matrices from

$$
\begin{bmatrix} \rho_w \\ \rho_v \end{bmatrix} = \begin{bmatrix} \hat{T} \hat{\chi}_s^+ \\ \mathcal{Y}_{p|p} \end{bmatrix} - \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} \hat{T} \hat{\chi}_s^+ \\ \mathcal{U}_{p|p} \end{bmatrix}. \quad (15)
$$

Finally, the noise covariance estimates are obtained from

$$
\begin{bmatrix} \hat{Q} & \hat{S} \\ \hat{S}^T & \hat{R} \end{bmatrix} = \begin{bmatrix} \rho_w \\ \rho_v \end{bmatrix} \begin{bmatrix} \rho_w^T & \rho_v^T \end{bmatrix}, \qquad (16)
$$

and can be used in the Kalman filter for state estimation in (4)–(5). It should
be noted that the identification of (Q, R, S) from data is not unique, since dif-
ferent noise covariance matrices can yield the same Kalman predictor gain and
thus state estimates [2]. In particular, the obtained estimates (Q̂, R̂, Ŝ) here are
consistent estimates of the noise covariances corresponding to the system in inno-
vation form and not directly of (Q, R, S); however, they yield the same Kalman
predictor gain and state estimates. The detailed proof of the latter is omitted
here for brevity, and can be found in [4].
## 4 Numerical validation

In this section the subspace-based approach for noise covariance estimation outlined in Section 3 is compared with the ALS method [1]. In what follows, the estimate of the noise covariance triplet obtained with the ALS approach is denoted ($\hat{Q}^{\text{ALS}}$, $\hat{R}^{\text{ALS}}$, $\hat{S}^{\text{ALS}}$). Both noise covariance estimation methods are validated on an example of a 6 DOF mechanical chain-like system. For any consistent set of units, the system is modeled with spring stiffness $k_1 = k_3 = k_5 = 100$ and $k_2 = k_4 = k_6 = 200$, mass $m_i = 1/20$ and a proportional damping matrix such that all modes have a damping ratio of 2%. The chain is illustrated in Figure 1. The system is excited by a white noise signal acting at all DOFs. An additional white noise with 1% of the standard deviation of the measured input is added to the input to emulate the process noise. The acceleration outputs are simulated at DOFs 1, 3 and 5 at a sampling frequency of 50 Hz, and white measurement noise with 5% of the standard deviation of the output is added to each response measurement. The corresponding ‘true’ noise covariance of the system are denoted by $(Q, R, S)$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. 6 DOF chain system sketch.</figcaption></figure>

A data set {$y_k$} at the output DOFs 1, 3 and 5 is simulated, where the true states {$x_k$} corresponding to the displacements at DOFs 2, 4 and 6 will be used as reference for state estimation. The data set {$y_k$} is used together with the inputs {$u_k$} from DOFs 3 and 4 to obtain estimates ($\hat{Q}$, $\hat{R}$, $\hat{S}$) with the method above.

The state prediction errors for the j-th state (normalized RMSE)

$$ \sqrt{\frac{\sum_{k} (\hat{x}_{k|k-1}^j - x_k^j)^2}{\sum_{k} (x_k^j)^2}} \quad (17) $$

are evaluated for the displacements at DOFs 2, 4 and 6 and averaged. The state prediction is based on the dataset {$y_k$} and the inputs {$u_k$} from above, the model matrices $(A, B, C, D)$ and the noise covariances

- true $(Q, R, S)$

- estimated with the subspace approach ($\hat{Q}$, $\hat{R}$, $\hat{S}$)

- estimated with ALS ($\hat{Q}^{\text{ALS}}$, $\hat{R}^{\text{ALS}}$, $\hat{S}^{\text{ALS}}$)

- a set of sampled $(Q_i, R_i, S_i)$ with a latin hypercube (LHC) [6]
<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison between the true and estimated states.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. LHS centered at true (Q, R, S) (left) and estimates ($\hat{Q}$, $\hat{R}$, $\hat{S}$) (right). Large latin hypercube range.</figcaption></figure>

For the LHC sampling, $i = 1, ..., 10,000$ samples are generated and centered either on the true (Q, R, S) or the subspace estimates ($\hat{Q}$, $\hat{R}$, $\hat{S}$). To ensure positive semidefiniteness of the sampled covariance, only a lower triangular matrix $L$ is generated so that $LL^T$ corresponds to the block matrix containing Q, R, S.

Figure 2 illustrates a comparison between the true states of the system, the states corresponding to the Kalman predictor obtained with the estimated noise covariance and the states corresponding to the Kalman predictor obtained with the LHC samples yielding the best and the worst RMSE. Both the left and the right parts of the figure show that the difference between the true states and the states related to the estimates of the noise covariance is small. Figure 3 illustrates the normalized RMSE from one data set for the estimated and true noise covariance, as well as RMSE histograms where the LHS is centered on the true noise covariance and on the subspace estimates. It can be seen that the theoretical noise covariance yields the smallest RMSE, which is in the range of several LHC samples. The subspace estimates of the noise covariance yield a smaller RMSE than the covariance obtained with the ALS approach and are within the 5% quantile of the distribution of RMSE of LHC samples.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Monte Carlo histograms of RMSE based on true, subspace-based and ALS noise covariance.</figcaption></figure>

The consistency of the subspace and ALS estimates is verified in a Monte Carlo experiment with 1000 samples. For each sample a new set of input/output data is generated, and a new set of {$Q, R, S$} triplet is estimated with the subspace and the ALS approaches. The resulting histograms of RMSE are compared with the histogram of RMSE stemming from the state predictions with the true noise covariance, which is illustrated in Figure 4.

It can be viewed that the subspace-based method on average yields a smaller RMSE than the ALS approach and, as expected, both yield errors that are significantly larger than RMSE resulting from using the true noise covariance.

# 5 Conclusions

In this work we reviewed a subspace-based approach for estimation of the noise covariance for Kalman filters and compared it to the state of art ALS approach in the virtual sensing context. Results show that the subspace-based method yields ($Q, R, S$) triplets which on average lead to more accurate state estimates than the ALS method. Further details and applications can be found in [4]. Future work will comprise statistical analysis of both approaches and an application of both to virtual sensing of a large-scale mechanical system.

# References

1. Åkesson, B.M., Jørgensen, J.B., Poulsen, N.K., Jørgensen, S.B.: A generalized autocovariance least-squares method for kalman filter tuning. Journal of Process Control **18**(7), 769–779 (2008)

2. Arnold, T.J., Rawlings, J.B.: Uniqueness conditions for ALS problems. IFAC-PapersOnLine **51**(20), 469–474 (2018), 6th IFAC Conference on Nonlinear Model Predictive Control NMPC 2018
3. Bavdekar, V.A., Deshpande, A.P., Patwardhan, S.C.: Identification of process and measurement noise covariance for state and parameter estimation using extended Kalman filter. Journal of Process Control **21**(4), 585–601 (2011)

4. Greś, S., Döhler, M., Dertimanis, V., Chatzi, E.: Subspace-based noise covariance estimation for Kalman filter in virtual sensing applications. Mechanical Systems and Signal Processing, submitted

5. Matisko, P., Havlena, V.: Noise covariance estimation for Kalman filter tuning using Bayesian approach and Monte Carlo. International Journal of Adaptive Control and Signal Processing **27**(11), 957–973 (2013)

6. McKay, M.D., Beckman, R.J., Conover, W.J.: A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics **42**(1), 55–61 (2000)

7. Mussot, V., Mercère, G., Dairay, T., Arvis, V., Vayssettes, J.: Noise covariance matrix estimation with subspace model identification for Kalman filtering. International Journal of Adaptive Control and Signal Processing **35**(4), 591–611 (2021)

8. Myers, K., Tapley, B.: Adaptive sequential estimation with unknown noise statistics. IEEE Transactions on Automatic Control **21**(4), 520–523 (1976)

9. Odelson, B.J., Rajamani, M.R., Rawlings, J.B.: A new autocovariance least-squares method for estimating noise covariances. Automatica **42**(2), 303–308 (2006)

10. van Overschee, P., de Moor, B.: Subspace Identification for Linear Systems. Springer, 1st edn. (1996)

11. Partovibakhsh, M., Liu, G.: An adaptive unscented Kalman filtering approach for online estimation of model parameters and state-of-charge of lithium-ion batteries for autonomous mobile robots. IEEE Transactions on Control Systems Technology **23**(1), 357–363 (2015)

12. Peeters, B., de Roeck, G.: Reference-based stochastic subspace identification for output-only modal analysis. Mechanical Systems and Signal Processing **13**(6), 855–878 (1999)

13. Reynders, E., de Roeck, G.: Reference-based combined deterministic-stochastic subspace identification for experimental and operational modal analysis. Mechanical Systems and Signal Processing **22**(3), 617–637 (2008)

14. Yuen, K.V., Liu, Y.S., Yan, W.J.: Estimation of time-varying noise parameters for unscented Kalman filter. Mechanical Systems and Signal Processing **180**, 109439 (2022)

15. Zagrobelny, M.A., Rawlings, J.B.: Identifying the uncertainty structure using maximum likelihood estimation. In: 2015 American Control Conference (ACC). pp. 422–427 (2015)