A Gaussian process latent force model for joint input-state estimation in linear structural systems

TL;DR

This paper introduces a Gaussian process latent force model for joint input-state estimation in linear structures, combining physical models with data-driven Gaussian processes, resulting in improved robustness and accuracy over traditional methods.

Contribution

The paper develops a novel GPLFM approach that integrates Gaussian processes with differential equations for enhanced joint estimation in structural systems.

Findings

Superior estimation accuracy compared to conventional Kalman filter methods

Robustness and numerical stability demonstrated in simulations

Hyperparameters estimated automatically via maximum likelihood

Abstract

The problem of combined state and input estimation of linear structural systems based on measured responses and a priori knowledge of structural model is considered. A novel methodology using Gaussian process latent force models is proposed to tackle the problem in a stochastic setting. Gaussian process latent force models (GPLFMs) are hybrid models that combine differential equations representing a physical system with data-driven non-parametric Gaussian process models. In this work, the unknown input forces acting on a structure are modelled as Gaussian processes with some chosen covariance functions which are combined with the mechanistic differential equation representing the structure to construct a GPLFM. The GPLFM is then conveniently formulated as an augmented stochastic state-space model with additional states representing the latent force components, and the joint input and state inference of the resulting model is implemented using Kalman filter. The augmented state-space model of GPLFM is shown as a generalization of the class of input-augmented state-space models, is proven observable, and is robust compared to conventional augmented formulations in terms of numerical stability. The hyperparameters governing the covariance functions are estimated using maximum likelihood optimization based on the observed data, thus overcoming the need for manual tuning of the hyperparameters by trial-and-error. To assess the performance of the proposed GPLFM method, several cases of state and input estimation are demonstrated using numerical simulations on a 10-dof shear building and a 76-storey ASCE benchmark office tower. Results obtained indicate the superior performance of the proposed approach over conventional Kalman filter based approaches.

Tables

Table 1. Table 1: Undamped natural frequencies and damping ratios of the 10-dof shear building

Mode	1	2	3	4	5	6	7	8	9	10
Frequency ($\mathrm{Hz}$)	1.19	3.54	5.81	7.96	9.92	11.67	13.15	14.34	15.21	15.74
Damping ratio (%)	0.86	0.78	1.05	1.35	1.64	1.90	2.13	2.31	2.44	2.52

Equations

$$\mathbf{M}\ddot{\mathbf{u}}(t)+\mathbf{C}\dot{\mathbf{u}}(t)+\mathbf{K}\mathbf{u}(t)={\mathbf{S}}_{p}p(t)-\mathbf{M}{\mathbf{S}}_g{\ddot{\mathbf{u}}}^g(t),$$

$$x(t)=\left[\begin{array}{c}\mathbf{u}(t)\\ \dot{\mathbf{u}}(t)\end{array}\right],f(t)=\left[\begin{array}{c}p(t)\\ {\ddot{\mathbf{u}}}^g(t)\end{array}\right],$$

$$\dot{x}(t)={\mathbf{A}}_{\mathrm{c}}x(t)+{\mathbf{B}}_{\mathrm{c}}f(t)$$

$${\mathbf{A}}_{\mathrm{c}}=\left[\begin{array}{cc}0& I\\ -{\mathbf{M}}^{-1}\mathbf{K}& -{\mathbf{M}}^{-1}\mathbf{C}\end{array}\right],{\mathbf{B}}_{\mathrm{c}}=\left[\begin{array}{cc}0& 0\\ {\mathbf{M}}^{-1}{\mathbf{S}}_p& -{\mathbf{S}}_g\end{array}\right],$$

$$y(t)=\left[\begin{array}{ccc}{\mathbf{S}}_{dis}& 0& 0\\ 0& {\mathbf{S}}_{vel}& 0\\ 0& 0& {\mathbf{S}}_{acc}\end{array}\right]\left[\begin{array}{c}\mathbf{u}(t)\\ \dot{\mathbf{u}}(t)\\ \ddot{\mathbf{u}}(t)\end{array}\right],$$

$$y(t)={\mathbf{G}}_{\mathrm{c}}x(t)+{\mathbf{J}}_{\mathrm{c}}f(t),$$

$${\mathbf{G}}_{\mathrm{c}}=\left[\begin{array}{cc}{\mathbf{S}}_{dis}& 0\\ 0& {\mathbf{S}}_{vel}\\ -{\mathbf{S}}_{acc}{\mathbf{M}}^{-1}\mathbf{K}& -{\mathbf{S}}_{acc}{\mathbf{M}}^{-1}\mathbf{C}\end{array}\right],{\mathbf{J}}_{\mathrm{c}}=\left[\begin{array}{cc}0& 0\\ 0& 0\\ {\mathbf{S}}_{acc}{\mathbf{M}}^{-1}{\mathbf{S}}_p& 0\end{array}\right].$$

$$y(t_k)={\mathbf{G}}_{c}x(t_k)+{\mathbf{J}}_{c}f(t_k)+v_k,$$

$$\begin{array}{rl}\dot{x}(t)& ={\mathbf{A}}_{\mathrm{c}}x(t)+{\mathbf{B}}_{\mathrm{c}}f(t),\\ y(t_k)& ={\mathbf{G}}_{c}x(t_k)+{\mathbf{J}}_{c}f(t_k)+v_k,\end{array}$$

$$\begin{array}{rl}{x}_k& =\mathbf{A}{x}_{k-1}+\mathbf{B}{f}_{k-1}+w_{k-1},\\ y_k& =\mathbf{G}{x}_k+\mathbf{J}{f}_k+v_k,\end{array}$$

$$\mathbb{E}\left[\left\{\begin{array}{c}{w}_k\\ v_k\end{array}\right\}\left\{\begin{array}{cc}{w}_l^T& v_l^T\end{array}\right\}\right]=\left[\begin{array}{cc}{Q}^x& 0\\ 0^T& R\end{array}\right]{\delta }_{kl},$$

$$g(\mathbf{\ell })\sim \mathrm{GP}\left(\mu (\mathbf{\ell }),\kappa (\mathbf{\ell },{\mathbf{\ell }}^{\prime };\theta )\right),$$

$$\begin{array}{rl}\mu (\mathbf{\ell })& =\mathbb{E}\left[g(\mathbf{\ell })\right],\\ \kappa (\mathbf{\ell },{\mathbf{\ell }}^{\prime };\theta )& =\mathbb{E}\left[\left(g(\mathbf{\ell })-\mu (\mathbf{\ell })\right)\left(g({\mathbf{\ell }}^{\prime })-\mu ({\mathbf{\ell }}^{\prime })\right)\right],\end{array}$$

$$\kappa (\mathbf{\ell },{\mathbf{\ell }}^{\prime })={\sigma }_g^2\exp \left[-\sum _{k=1}^d\frac{{\left(\mathbf{\ell }(k)-{\mathbf{\ell }}^{\prime }(k)\right)}^2}{2l_k^2}\right],$$

$$y=g(\mathbf{\ell })+\nu ,$$

$$y_k=g(t_k)+{\nu }_k,{\nu }_k\sim \mathcal{N}\left(0,{\sigma }_n^2\right).$$

$$g(t)\sim \mathrm{GP}(0,k(t,t^{\prime };\theta )),$$

$$\left[\begin{array}{c}y\\ g(t_{\ast })\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}\mathbf{K}(\mathbf{t},\mathbf{t};\theta )+{\sigma }_n^{2}I& \mathbf{k}(\mathbf{t},t_{\ast };\theta )\\ {\mathbf{k}}^T(\mathbf{t},t_{\ast };\theta )& k(t,t_{\ast };\theta )\end{array}\right]\right),$$

$$p(g(t_{\ast })|y,\theta )=\mathcal{N}\left(g(t_{\ast })|{\mu }_{\mathrm{GP}}(t_{\ast }),{\sigma }_{\mathrm{GP}}^2(t_{\ast })\right),$$

$$\begin{array}{rl}{\mu }_{\mathrm{GP}}(t_{\ast })& ={\mathbf{k}}^T(\mathbf{t},t_{\ast };\theta ){\left[\mathbf{K}(\mathbf{t},\mathbf{t};\theta )+{\sigma }_n^{2}I\right]}^{-1}y,\\ {\sigma }_{\mathrm{GP}}^2(t_{\ast })& =k(t_{\ast },t_{\ast };\theta )-{\mathbf{k}}^T(\mathbf{t},t_{\ast };\theta ){\left[\mathbf{K}(\mathbf{t},\mathbf{t};\theta )+{\sigma }_n^{2}I\right]}^{-1}\mathbf{k}(\mathbf{t},t_{\ast };\theta ).\end{array}$$

$$\begin{array}{rl}\dot{x}(t)& ={\mathbf{A}}_{\mathrm{c}}x(t)+{\mathbf{B}}_{\mathrm{c}}f(t),\\ y(t_k)& ={\mathbf{G}}_{c}x(t_k)+{\mathbf{J}}_{c}f(t_k)+v_k,\end{array}$$

$$f^{(j)}(t)\sim \mathrm{GP}(0,{\kappa }_j(t,t^{\prime })),$$

$$x(t)=\mathbf{\Psi }(t)x_0+{\int }_0^t\mathbf{\Psi }(t-s){\mathbf{B}}_{\mathrm{c}}f(s)ds,$$

$${\mathbf{K}}_{xx}(t,t^{\prime })=\mathbb{E}\left[x(t)x(t^{\prime }{)}^T\right]=\mathbf{\Psi }(t)P_0^x\mathbf{\Psi }(t{)}^T+{\int }_0^{t^{\prime }}{\int }_0^t\mathbf{\Psi }(t-s){\mathbf{B}}_{\mathrm{c}}{\mathbf{K}}_{ff}(s,s^{\prime }){{\mathbf{B}}_{\mathrm{c}}}^T\mathbf{\Psi }(t^{\prime }-s^{\prime }{)}^{T}dsd{s}^{\prime },$$

$${\mathbf{K}}_{ff}(t,t^{\prime })=\mathbb{E}\left[f(t)f(t^{\prime }{)}^T\right]=\text{diag}\left[{\kappa }_1(t,t^{\prime }),\dots ,{\kappa }_{n_f}(t,t^{\prime })\right].$$

$$x(t)\sim \mathrm{GP}(0,{\mathbf{K}}_{xx}(t,t^{\prime })).$$

$$y_k=g(t_k)+v_k$$

$$g(t_k)={\mathbf{G}}_{c}x(t_k)+{\mathbf{J}}_{c}f(t_k)=\underset{{\mathbf{G}}_a}{\underbrace{\left[\begin{array}{cc}{\mathbf{G}}_c& {\mathbf{J}}_c\end{array}\right]}}\underset{\eta (t)}{\underbrace{\left[\begin{array}{c}x(t)\\ f(t)\end{array}\right]}}={\mathbf{G}}_a\eta (t).$$

$$g(t)\sim \mathrm{GP}(0,{\mathbf{K}}_{gg}(t,t^{\prime })),$$

$${\mathbf{K}}_{gg}(t,t^{\prime })={\mathbf{G}}_a{\mathbf{K}}_{\eta \eta }(t,t^{\prime }){\mathbf{G}}_a^T,$$

Full text

K****eywords Input estimation  $\cdot$ state estimation  $\cdot$ force identification  $\cdot$ latent force models  $\cdot$ Gaussian process  $\cdot$ linear system  $\cdot$ time-invariant

1 Introduction

Ordinary differential equations relating the applied input forces to the system state variables are commonly employed in modeling structural dynamic systems. Direct measurements of the system state variables are not always possible; instead, noisy observations as a function of the system state variables are obtained at discrete time instants. State estimation addresses the problem of obtaining optimal estimates of the unobserved system states from noisy measurements given the knowledge of the system model and the dynamic input forces. This problem is at the core of many engineering applications ranging from condition monitoring and prediction of stresses, to feedback control and performance evaluation. The Kalman filter and its variants have been widely used for state estimation of a broad class of both linear and nonlinear dynamical systems [1, 2, 3]. A significant challenge in dynamical state estimation arises due to inadequate knowledge of the inputs acting on the system. In many civil engineering systems such as buildings or bridges, measuring the operational external forces (wind loads, vehicular loads, etc.) is not practical. A commonly adopted approach in such cases of unmeasured input forces is to assume that the input force is a zero mean white noise process and then apply filtering techniques for state estimation; however, in many cases the actual inputs may violate this white noise assumption and could potentially yield poor estimation results. Thus, there is a need to inversely determine the unmeasured input forces. Input estimation, also known as force identification, involves reconstructing the latent or unmeasured dynamic input forces given a limited number of sensor measurements and knowledge of the system model —obtained either from a finite element model or inferred a priori using system identification methods. Cases where input estimation and state estimation are performed together are referred to as joint input-state estimation.

Originally, force identification problems were treated separately from state estimation problems, and the earliest algorithms used frequency response functions to inversely compute the unknown forces [4, 5]. Subsequently, a wide variety of force identification methods were proposed, both in time-domain and frequency-domain (see [6] for a recent survey on force reconstruction techniques). Among time-domain approaches, two distinct modelling frameworks exist: deterministic [7, 8] and stochastic [9, 10, 11]. The stochastic time-domain approaches differ from their deterministic counterparts in the manner that uncertainty in the form of noise is incorporated in the system and measurement model, and the state and/or input estimates so obtained under stochastic framework are characterized by probability distributions.

Recently, there has been a considerable amount of interest in developing and applying stochastic joint input-state estimation algorithms. Gillijns and De Moor [12] proposed a minimum-variance unbiased filter for joint input and state estimation of linear systems without a direct transmission term, where the input estimation and state estimation are performed in two independent, sequential steps. This algorithm was extended to include the direct transmission term [13] which was later applied to structural dynamics by Lourens et al. [14] for use with reduced-order models (i.e., models obtained using a relatively small number of modes). Lourens et al. [15] also proposed an augmented Kalman filter (AKF) that appends the unknown forces to the state vector to jointly estimate the dynamic forces and system states using a classical Kalman filter. This method results in spurious low frequency components in the force and displacement estimates when only acceleration measurements are used. An analytical investigation performed in [16] showed that the augmented state-space formulation of AKF is susceptible to numerical instability and un-observability issues when only acceleration measurements are used. In an effort to overcome these shortcomings, Naets et al. [16] proposed the use of artificial displacement measurements to stabilize systems with only acceleration measurements. This algorithm, known as AKF with dummy measurements (AKFdm), is similar to the one proposed by Chatzi and Fuggini [17] for alleviating the low frequency drifts in displacement estimates needed for structural monitoring purposes. Azam et al. [18] proposed a dual Kalman filter (DKF) for state and input estimation using only acceleration measurements. Unlike the AKF formulation where the input and states are jointly predicted and updated, the DKF follows a successive structure of implementation, in that, the input prediction and update are followed by the state prediction and update. Due to the separation of the estimation process in two successive steps, the DKF typically overcomes the issue of un-observability. Maes et al. [19] proposed a modified AKF which directly accounts for the correlation between the process noise and the measurement noise.

An attractive feature of the above Kalman filtering-based joint input-state estimation methods lies in their amenability for online implementation. That said, their accuracy in force reconstruction critically depends on prior tuning of the covariance matrices associated with the unknown inputs. For example, in the AKF and DKF algorithms, the covariance matrix for the unknown input is tuned while in the AKFdm, an additional covariance matrix for the dummy displacements must be tuned in addition to the input covariance. In the absence of expert knowledge of the covariance matrices, a trial-and-error approach or a heuristic rule of thumb is often adopted to tune the input covariance matrices prior to implementation. This may not be possible when there is very little prior knowledge about the input and/or the states of the system. The tacit approach followed in most cases is to have an offline tuning phase —where some response data collected over a certain time duration is used to obtain optimal values of the covariance matrices —prior to the online estimation phase with Kalman filtering.

In 2009, Alvarez et al. [20] proposed a Bayesian learning technique termed latent force model (LFM), that allows flexible modelling of unknown inputs (or latent forces) in a system. The unknown or latent forces of a system are modelled using non-parametric Gaussian process (GP) model [21] with some chosen covariance functions. The choice of covariance function provides flexibility to incorporate generic prior knowledge about the behavior of the inputs (e.g. random, periodic, smoothly varying). A Gaussian process based LFM, GPLFM in short, is then constructed incorporating the physics-based mathematical model of the system (e.g. differential equations of the system) within the GP covariance functions to infer inputs to the system. The behavior of the latent forces depend on the parameters (also called hyperparameters) of the chosen GP covariance functions which constitute the tunable parameters in this model. These hyperparameters can be set based on expert knowledge, or they can be tuned by optimization using observed data. Hartikainen and Särkkä [22] and more recently Särkkä et al. [23] advocated a state-space approach to GPLFM inference wherein all information about the system and the GP force model is fully captured in an augmented state-space representation, thus making the inference amenable to use with Bayesian filtering and smoothing methods. Once the hyperparameters of GPLFM are computed, a GPLFM can be adapted for online implementation with Kalman filter.

This paper introduces GPLFM for joint input and state estimation of linear time-invariant structural systems. To the authors’ knowledge, this is the first application of GPLFM in the domain of structural dynamics. Additionally, the following are the key contributions of this paper:

The GPLFM is shown to be a generalization of the class of input-augmented state-space models underlying some popular Kalman filter-based joint input-state estimation approaches [15, 16, 19]. The augmented state-space model used in AKF [15] is obtained as a special case of GPLFM. 2. 2.

Unlike the AKF, the GPLFM formulation is analytically proved to be observable and is shown to be robust against drift in force estimation even when used with only acceleration measurements. 3. 3.

A maximum likelihood optimization based on the observed data is used to compute the optimal hyperparameters of the GPLFM, which eliminates the need for tuning the parameters heuristically or through trial-and-error.

An in-depth numerical study is also presented comparing the performance of the proposed GPLFM algorithm with the existing AKF, AKFdm, and DKF algorithms.

The paper is organized as follows. Section 2 presents the mathematical model of the structure and outlines the objective of this study. Next, a concise background on GPLFM and its use in regression is provided in Section 3. In Section 4, the GPLFM is formulated into a linear state-space model, and a procedure for joint posterior inference on the inputs and states is discussed. Section 4.5 presents GPLFM as a generalization of the augmented state-space models used in joint input-state estimation and provides proofs of observability and stability of the GPLFM. Section 6 encompasses numerical studies using a 10-dof shear building, illustrating the performance of GPLFM under different excitation and measurement scenarios. The paper is finally concluded with an application of the proposed method on a 76-storey ASCE benchmark office tower for joint input and state estimation.

2 Problem Formulation

2.1 Mathematical model of structural system

The equation of motion of a $\operatorname{\mathit{n}_{\text{dof}}}$ degrees-of-freedom (dofs) linear structural system under forced excitation including ground motions can be represented by the following second order ordinary differential equation (obtained after discretization in space):

[TABLE]

where, $\bm{{u}}(t)\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\text{dof}}}}$ is the vector of displacements at the degrees of freedom and $\mathbf{{M}}$, $\mathbf{{C}}$ and $\mathbf{{K}}\in\mathbb{R}^{\operatorname{\mathit{n}_{\text{dof}}}\times\operatorname{\mathit{n}_{\text{dof}}}}$ represent the mass, damping, and stiffness matrices of the structural system, respectively. The external forces acting on the structure are represented by a combination of the external loads acting on the structure $\mathbf{{S}}_{p}\operatorname{\bm{{p}}}(t)$ and the forces generated due to earthquake ground motion $\mathbf{{M}}\mathbf{{S}}_{g}\operatorname{\ddot{\bm{{u}}}}^{g}(t)$. The first term denotes the product of a load influence matrix $\mathbf{{S}}_{p}\in\mathbb{R}^{\operatorname{\mathit{n}_{\text{dof}}}\times\operatorname{\mathit{n}_{\mathit{p}}}}$ and the vector $\operatorname{\bm{{p}}}(t)\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{p}}}}$ representing the $\operatorname{\mathit{n}_{\mathit{p}}}$ external load time histories while the second term represents the forces generated in the structure due to earthquake ground motion $\bm{{u}}^{g}(t)\in\operatorname{\mathbb{R}}^{n_{g}}$ where $\mathbf{{S}}_{g}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\text{dof}}}\times n_{g}}$ is the matrix of influence coefficients which specifies the dofs affected by the ground motion. In cases where the size of the structural model is very large, it is common to use modal reduced-order models.

Defining a state vector $\operatorname{\bm{{x}}}(t)$ and a force vector $\operatorname{\bm{{f}}}(t)$ as

[TABLE]

where $\operatorname{\bm{{x}}}(t)\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{x}}}}$ and $\operatorname{\bm{{f}}}(t)\in\operatorname{\mathbb{R}}^{n_{f}}$; $\operatorname{\mathit{n}_{\mathit{x}}}=2\operatorname{\mathit{n}_{\text{dof}}}$ and $n_{f}=n_{p}+n_{g}$, the set of second-order differential equations in Equation 1 can be converted to a set of first-order differential equations called as the continuous-time state-space equation

[TABLE]

where the continuous-time system matrices $\operatorname{\mathbf{{A}}_{c}}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{x}}}\times\operatorname{\mathit{n}_{\mathit{x}}}}$ and $\operatorname{\mathbf{{B}}_{c}}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{x}}}\times\operatorname{\mathit{n}_{\mathit{f}}}}$ are defined as

[TABLE]

where $\operatorname{\mathbf{{I}}}$ is the identity matrix of appropriate dimension. Consider next the measurement equation, where it is assumed that a combination of displacements, velocities and accelerations can be measured. Hence the output vector $\operatorname{\bm{{y}}}(t)$, containing $\operatorname{\mathit{n}_{\mathit{y}}}$ measured quantities, assumes the following form

[TABLE]

where, $\mathbf{{S}}_{dis}$, $\mathbf{{S}}_{vel}$ and $\mathbf{{S}}_{acc}$ are the selection matrices for displacements, velocities and accelerations, respectively. Using Equation 1 and the defined state and force vectors, Equation 5 can be written in the state-space form as

[TABLE]

where, the output influence matrix $\operatorname{\mathbf{{G}}_{c}}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{y}}}\times\operatorname{\mathit{n}_{\mathit{x}}}}$ and direct transmission matrix $\operatorname{\mathbf{{J}}_{c}}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{y}}}\times\operatorname{\mathit{n}_{\mathit{f}}}}$ are

[TABLE]

Note that the second column of the direct transmission term $\operatorname{\mathbf{{J}}_{c}}$ corresponds to the earthquake ground motion input and takes zero values as the direct contribution of earthquake ground motion vanishes when measuring absolute accelerations (see derivation in D)

Equations 3 and 6 form the continuous-time state-space model (SSM) for the system described by Equation 1. In practical applications, continuous time outputs $\operatorname{\bm{{y}}}(t)$ are of course not observed, instead noisy measurements are obtained via sampling of system responses at discrete time instants. Assuming a sampling interval of $\Delta t$, discrete time instances are defined at $t_{k}=k\Delta t$ for $k=1,\ldots,N$, and a discrete measurement model can be expressed as

[TABLE]

where $\operatorname{\bm{{v}}}_{k}$ is the measurement noise vector. Combined together, Equations 3 and 8 represent a SSM with continuous-time dynamics and discrete-time measurements known as the continuous-discrete state-space model (see [24], pg170 or [25], pg93)

[TABLE]

For numerical implementation, the continuous-time state-space form in Equations 3 and 6 is converted to the discrete-time state-space following the zero-order-hold assumption

[TABLE]

where $\mathbf{{A}}=\exp(\operatorname{\mathbf{{A}}_{c}}\Delta t)$, $\mathbf{{B}}=[\mathbf{{A}}-\operatorname{\mathbf{{I}}}]\operatorname{\mathbf{{A}}_{c}}^{-1}\operatorname{\mathbf{{B}}_{c}}$, $\mathbf{{G}}=\operatorname{\mathbf{{G}}_{c}}$, $\mathbf{{J}}=\operatorname{\mathbf{{J}}_{c}}$. The terms $\operatorname{\bm{{w}}}_{k-1}$ and $\operatorname{\bm{{v}}}_{k}$ are the process and the measurement noise vectors, typically added to to account for modelling errors and measurement errors, respectively. They are assumed zero-mean Gaussian white noise with covariances

[TABLE]

with $\operatorname{\mathbf{{R}}}\succ\operatorname{\mathbf{{0}}}$ and $\operatorname{\mathbf{{Q}}}^{x}\succcurlyeq\operatorname{\mathbf{{0}}}$.

2.2 Objective

Having elaborated on the mathematical model of structure of interest, attention is now focused on the problem of stochastic joint input and state estimation. Here, both the inputs $\operatorname{\bm{{f}}}_{k}$ and the system states $\operatorname{\bm{{x}}}_{k}$ are taken to be unknown sequence of Gaussian random variables, and the goal lies in inferring the joint posterior distribution of the inputs and states from a set of noisy measurements $\operatorname{\bm{{y}}}_{k}$ given the a priori knowledge of the structural model parameters. In pursuit of this goal, the objective of this paper is to implement a Gaussian process latent force model which improves upon existing methods by reducing the dependency on manual tuning of covariance matrices associated with the unknown inputs and also provides numerical stability with respect to observability of the augmented state-space formulation when used with only acceleration measurements. The location of the inputs is assumed to be known, however studies on input localization can be found in [26, 27, 28].

3 Background on Gaussian process latent force models

3.1 Brief overview of Gaussian processes

Gaussian processes (GPs) [29, 21, 30] are a class of stochastic processes that provide a paradigm for specifying probability distributions over functions. Gaussian processes have been widely studied and used in many different fields such as signal processing [31], geostatistics [32, 33] and inverse problems [30]. For instance, in geostatistics literature Gaussian process regression is known as Kriging. A Gaussian process is a generalization of the Gaussian probability distribution. While a probability distribution is defined over random variables which are scalars or vectors (for multivariate distributions), a stochastic process governs the properties of functions. Consider an independent variable $\bm{{\ell}}\in\operatorname{\mathbb{R}}^{d}$ and a function $g(\bm{{\ell}})$ such that $g:\operatorname{\mathbb{R}}^{d}\rightarrow\operatorname{\mathbb{R}}$. Then a GP defined over $g(\bm{{\ell}})$ with the mean $\mu(\bm{{\ell}})$ and covariance function $\kappa(\bm{{\ell}},\bm{{\ell}}^{\prime};\operatorname{\bm{{\theta}}})$ is as follows

[TABLE]

[TABLE]

where $\operatorname{\bm{{\theta}}}$ denotes the hyperparameters of the covariance function $\kappa$. The choice of the covariance function $\kappa$ allows encoding any prior knowledge about $g(\bm{{\ell}})$ (e.g., periodicity, linearity, smoothness), and can accommodate approximation of arbitrarily complex functions [21]. The notation in Equation 12 implies that any finite collection of function values has a joint multivariate Gaussian distribution, that is $\left(g(\bm{{\ell}}_{1}),g(\bm{{\ell}}_{2}),\ldots,g(\bm{{\ell}}_{N})\right)\sim\mathcal{N}(\bm{{\mu}},\mathbf{{K}})$, where $\bm{{\mu}}=\left[\mu({\bm{{\ell}}_{1}}),\ldots,\mu({\bm{{\ell}}_{N}})\right]^{T}$ is the mean vector and $\mathbf{{K}}$ is the covariance matrix with $\mathbf{{K}}(i,j)=\kappa(\bm{{\ell}}_{i},\bm{{\ell}}_{j})$ for $i,j=1,2,\ldots,N$. If no prior information is available about the mean function, it is generally set to zero, i.e. $\mu(\bm{{\ell}})=0$. However, for the covariance function, any function $\kappa(\bm{{\ell}},\bm{{\ell}}^{\prime})$ that generates a positive, semi-definite, covariance matrix $\mathbf{{K}}$ is a valid covariance function. For example, the following squared exponential covariance function has the form

[TABLE]

where $\operatorname{\bm{{\theta}}}=\left\{\sigma_{g},l_{1},\ldots,l_{d}\right\}$ are the hyperparameters of the covariance function. The mean and covariance functions dictate how the random functions $g(\bm{{\ell}})$ behave on average and how the different points in the input space vary with respect to each other. The covariance function thus encodes a correlation structure which introduces inter-dependencies between function values at different inputs.

3.2 Gaussian process regression

GPs are often used to solve regression problems [34, 35, 36]. In general, regression problems are concerned with the estimation of values of a dependent variable $g(\bm{{\ell}})$ observed at certain values of an independent variable $\bm{{\ell}}$, given a set of noisy measurements $y$. The relation can be written as

[TABLE]

where $\nu$ is noise. In other words, regression involves estimating the latent (unobserved) function $g(\bm{{\ell}})$ that will enable prediction of $\hat{y}$ at new values of $\bm{{\ell}}$.

One perspective of looking at GPs representing the input-output relation is called the function-space view, given in Rasmussen and Williams [21]. Unlike traditional modelling approaches which rely on fitting a parameterized mathematical form to approximate the input-output relation $g(\bm{{\ell}})$, a GP does not assume any explicit form of $g(\bm{{\ell}})$, rather a prior belief (in the form of the mean function and the covariance function) is placed on the space of all possible functions $g(\bm{{\ell}})$. GPs belongs to the class of ‘non-parametric’ models because the number of parameters in the model is not fixed, but rather determined by the number of data points. Upon data collection, the posterior distribution over $g(\bm{{\ell}})$ is updated according to Bayes’ rule. All the candidate functions represented by GPs can be used to express linear or nonlinear relationships between the input $\bm{{\ell}}$ and the output $y$.

It is to be remarked at this point that the input dimension $d$ (dimension of the independent variable $\bm{{\ell}}$) is usually greater than 1, however in this work, the primary focus is on doing inference using time-series data where time is the only independent input variable. Thus $d=1$, and $\bm{{\ell}}$ is hereafter replaced by $t$ — referred to as the time domain. Examples where input dimension is greater that one ($d>1$) are situations where the input variables may represent time and one-dimensional space ($d=2$), or space in two dimensions ($d=2$), or even time with two-dimensional space ($d=3$).

Now consider the prediction of the value of $g(t)$ at a test point $t_{*}$, based on a previously measured set of outputs $\mathcal{D}=\left\{(t_{1},y_{1}),\ldots,(t_{N},y_{N})\right\}$ which are corrupted with white Gaussian noise:

[TABLE]

Assume that $g(t)$ is modelled as a Gaussian process with zero mean and a covariance function, that is

[TABLE]

where $k(t,t^{\prime};\operatorname{\bm{{\theta}}})$ is the covariance function with hyperparameters $\operatorname{\bm{{\theta}}}$. Under the Gaussian process assumption, the joint distribution between $y_{1},\ldots,y_{N}$ and $g(t_{*})$ is Gaussian, and can be expressed as

[TABLE]

where $\bm{{t}}=\left[t_{1},\ldots,t_{N}\right]^{T}$ is the vector consisting of time points, $\operatorname{\bm{{y}}}=\left[y_{1},\ldots,y_{N}\right]^{T}$ is a vector of measured outputs, $\mathbf{{K}}(\bm{{t}},\bm{{t}};\operatorname{\bm{{\theta}}})\in\operatorname{\mathbb{R}}^{N\times N}$ is the covariance matrix comprising elements $\left\{k(t_{i},t_{j};\operatorname{\bm{{\theta}}})\right\}_{i,j=1}^{N}$, and $\bm{{k}}(\bm{{t}},t_{*};\operatorname{\bm{{\theta}}})\in\operatorname{\mathbb{R}}^{n}$ is a vector with elements $\left\{k(t_{i},t_{*};\operatorname{\bm{{\theta}}})\right\}_{i=1}^{N}$.

Using the properties of multivariate Gaussian distributions, it can be shown that the posterior distribution of $g(t_{*})$ given the dataset $\operatorname{\bm{{y}}}$ and hyperparameters $\operatorname{\bm{{\theta}}}$ is also Gaussian

[TABLE]

with mean and variance

[TABLE]

3.3 Latent force models using Gaussian processes

Latent force models [20, 37] are a hybrid approach of modelling which combines mechanistic models (i.e. models based on physical laws) with non-parametric components such as GPs. Consider the following linear continuous-time state-space equations obtained from Equations 3 and 8

[TABLE]

with all components of the unknown input vector $\operatorname{\bm{{f}}}(t)$ being modelled as zero-mean independent Gaussian processes

[TABLE]

Here $f^{(j)}(t)$ is the $j$th component of $\operatorname{\bm{{f}}}(t)$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$. The assumption of zero mean leads to no loss of generality. Note that the hyperparameters $\operatorname{\bm{{\theta}}}$ of the covariance function have been omitted for notational compactness. The model, given by Equations 21 and 22, leads to a (linear) latent force model, in which the basic idea is to combine the mechanistic state-space model with unknown (latent) forcing functions $\operatorname{\bm{{f}}}(t)$ modelled by non-parametric Gaussian processes.

The solution of the state in the differential equation in Equation 21 can be obtained as

[TABLE]

where $\bm{\Psi}(\tau)=e^{\operatorname{\mathbf{{A}}_{c}}\tau}$ denotes the matrix exponential and $\operatorname{\bm{{x}}}_{0}$ is the initial condition of state at time $t=0$. The state covariance matrix function $\mathbf{{K}}_{xx}$ can then be calculated as

[TABLE]

where $\operatorname{\mathbf{{P}}}_{0}^{x}$ is the prior covariance matrix of $\operatorname{\bm{{x}}}(t)$ and $\mathbf{{K}}_{ff}(t,t^{\prime})$ is the joint covariance matrix of all latent forcing functions between time instants $t$ and $t^{\prime}$

[TABLE]

$\mathbf{{K}}_{ff}(t,t^{\prime})$ is a diagonal matrix due the assumption of independence across all force components. Observe that the states $\operatorname{\bm{{x}}}(t)$ form a multi-dimensional GP

[TABLE]

and thus both $\operatorname{\bm{{f}}}(t)$ and $\operatorname{\bm{{x}}}(t)$ are GPs. Next the measurement equation is rewritten as

[TABLE]

with

[TABLE]

Since the linear combination of Gaussian processes is a Gaussian process, $\operatorname{\bm{{g}}}(t)$ is a Gaussian process, and can be expressed as

[TABLE]

where

[TABLE]

and $\mathbf{{K}}_{ff}(t,t^{\prime})$ and $\mathbf{{K}}_{xx}(t,t^{\prime})$ are defined by Equations 25 and 24 respectively. A derivation for $\mathbf{{K}}_{xf}(t,t^{\prime})$ can be found in C. Thus, one can replace the problem defined by Equations 21, and 22 with an equivalent Gaussian process regression problem

[TABLE]

This approach, introduced by Alvarez et al. [20, 37] reduces the latent force model in Equation 21 and 22 to a GP regression problem as expressed by Equation 31. However, the approach requires analytical computation of the covariance matrix functions $\mathbf{{K}}_{xx}$, which may not always be obtained in a closed form and may need computation via numerical integration. In the next section, a formulation to convert covariance functions to state-space form for a certain class of GP covariance functions is considered. In that formulation, numerical integration is avoided altogether and replaced with matrix exponential computations and solutions to Lyapunov equations.

4 GPLFM as stochastic SSMs for sequential inference

This section discusses how temporal GPs with a certain class of covariance functions can be reformulated as SSMs and demonstrates how they can be combined with linear LFMs for sequential posterior inference of latent inputs and states. It must be emphasized at this point that only stationary covariance functions which have the form $\kappa(t,t^{\prime})=\kappa(\tau)$, where $\tau=t-t^{\prime}$, are considered in this work.

4.1 State space representation of temporal GPs

Consider a scalar-valued GP with stationary covariance function as follows

[TABLE]

It has been shown [38, 22, 39, 40] that $h(t)$ can be represented as the output of a scalar linear time-invariant system driven by white noise

[TABLE]

where $a_{0},\ldots,a_{m-1}$ are known constants and $w(t)$ is a white noise process with spectral density $S_{w}(\omega)=\sigma_{w}$. Equation 33 can be written in state-space form

[TABLE]

where the state $\operatorname{\bm{{z}}}(t)=\begin{pmatrix}h(t)&\frac{dh(t)}{dt}&\ldots&\frac{d^{m-1}h(t)}{dt^{m-1}}\end{pmatrix}^{T}$ and the matrices $\operatorname{\mathbf{{F}}}_{cf}\in\operatorname{\mathbb{R}}^{m\times m}$, $\mathbf{{L}}_{cf}\in\operatorname{\mathbb{R}}^{m\times 1}$ and $\mathbf{{H}}_{cf}\in\operatorname{\mathbb{R}}^{1\times m}$ are given by:

[TABLE]

The discrete-time form of the LTI model in Equation 34 is given by

[TABLE]

where $\operatorname{\mathbf{{F}}}=e^{\operatorname{\mathbf{{F}}}_{cf}\Delta t}$, $\mathbf{{L}}=(\operatorname{\mathbf{{F}}}-\operatorname{\mathbf{{I}}})\operatorname{\mathbf{{F}}}_{cf}^{-1}\mathbf{{L}}_{cf}$ and $\mathbf{{H}}=\mathbf{{H}}_{cf}$. This formulation allows GP regression problems to be solved sequentially with Kalman filtering and smoothing techniques. The iteration can be started with initial state $\operatorname{\bm{{z}}}_{0}\sim\mathcal{N}(\operatorname{\mathbf{{0}}},\operatorname{\mathbf{{P}}}^{cf}_{0})$. Note that the model (and the corresponding covariance function) is stationary, and therefore $\operatorname{\bm{{z}}}_{k}$ has a steady state solution distributed as $\operatorname{\bm{{z}}}_{k}\sim\mathcal{N}(\operatorname{\mathbf{{0}}},\operatorname{\mathbf{{P}}}^{cf}_{\infty})$, where the steady state covariance matrix $\operatorname{\mathbf{{P}}}^{cf}_{\infty}$ can be obtained as the solution of continuous-time Riccati solution

[TABLE]

Since steady state is invariant with respect to time, the initial covariance matrix is set equal to the steady state covariance, $\operatorname{\mathbf{{P}}}^{cf}_{0}\coloneqq\operatorname{\mathbf{{P}}}^{cf}_{\infty}$.

Using the SSM given by Equation 34, the power spectral density of $h(t)$ can be computed as

[TABLE]

The stationary covariance function for $h(t)$ is related to its spectral density through the inverse Fourier transform

[TABLE]

In terms of the state-space matrices, $\kappa(\tau)$ can be calculated using [41]

[TABLE]

where $\bm{\Phi}(\tau)=e^{\operatorname{\mathbf{{F}}}_{cf}\tau}$.

Hartikainen and Särkkä [38] showed that it is possible to form $\operatorname{\mathbf{{F}}}_{cf}$, $\mathbf{{L}}_{cf}$ and $S_{h}(\omega)$ such that $h(t)$ has the desired covariance function $\kappa(\tau)$ only if the spectral density of $h(t)$ has a proper rational form

[TABLE]

By applying spectral factorization, the spectral density $S_{h}(\omega)$ can be expressed similar to Equation 38

[TABLE]

where $H(i\omega)$ is a stable rational transfer function – has poles in the left half plane. The procedure to convert a covariance function $\kappa(\tau)$ into SSM is outlined as follows [39]:

• •

The spectral density $S_{h}(\omega)$ is computed by taking Fourier transform of $\kappa(\tau)$

• •

Spectral factorization of $S_{h}(\omega)$ is used to find a stable rational transfer function $H(i\omega)$ and $\sigma_{w}$ (as shown in Equation 42). If $S_{h}(\omega)$ is not rational as shown in Equation 41, then an approximation is used to represent it in a rational form.

• •

Finally, the transfer function model is converted into an equivalent SSM driven by a scalar-valued white noise process with spectral density $\sigma_{w}$

Many GPs having stationary covariance functions can be expressed into state-space form as in Equation 36 with model matrices defined by $\operatorname{\mathbf{{F}}}_{cf}$, $\mathbf{{L}}_{cf}$, $\mathbf{{H}}$, $\sigma_{w}$ and $\operatorname{\mathbf{{P}}}^{cf}_{0}$. It should be noted that non-stationary covariance functions can also be converted into state-space forms [42], however the Fourier transform approach is no longer applicable.

4.2 GP-LFM in joint state-space form

Consider once again the latent force model given by Equations 21 and 22, rewritten for convenience,

[TABLE]

where each component, $f^{(j)}(t)\in\operatorname{\mathbb{R}}$, of the latent force vector $\operatorname{\bm{{f}}}(t)$ is modelled as a GP. For each component, given a stationary covariance function $\kappa_{j}(\tau;\operatorname{\bm{{\theta}}}_{j})$ with hyperparameters $\operatorname{\bm{{\theta}}}_{j}$, a state-space representation can be constructed using the above procedure. The state-space representation for the $j$th component of $\operatorname{\bm{{f}}}(t)$, is given by

[TABLE]

where $\operatorname{\bm{{z}}}^{(j)}(t)\in\operatorname{\mathbb{R}}^{m_{j}}$ is the state vector for $j$th latent force, $w^{(j)}(t)$ is a zero mean Gaussian process with spectral density $\sigma_{w}^{(j)}$, and $\operatorname{\mathbf{{F}}}^{(j)}\in\operatorname{\mathbb{R}}^{m_{j}\times m_{j}}$, $\mathbf{{L}}^{(j)}\in\operatorname{\mathbb{R}}^{m_{j}\times 1}$, $\mathbf{{H}}^{(j)}\in\operatorname{\mathbb{R}}^{1\times m_{j}}$ are the state-space matrices with $m_{j}$ being the dimension of the state vector $\operatorname{\bm{{z}}}^{(j)}(t)$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$. Formally, the linear time-invariant system represented in Equation 44 is a stochastic differential equation where $w^{(j)}(t)$ represents a Brownian motion, and the mathematical description is given via stochastic integrals (see e.g. [24] for an account of this). The state-space form of the LFM can be obtained by combining the component GP SSMs with the system SSM to yield the following augmented model:

[TABLE]

where $\bm{{b}}_{1},\bm{{b}}_{2},\ldots,\bm{{b}}_{\operatorname{\mathit{n}_{\mathit{f}}}}$ are the columns of $\operatorname{\mathbf{{B}}_{c}}$ matrix and $\bm{{j}}_{1},\bm{{j}}_{2},\ldots,\bm{{j}}_{\operatorname{\mathit{n}_{\mathit{f}}}}$ are the columns of $\mathbf{{J}}_{c}$ matrix and vector $\tilde{\operatorname{\bm{{w}}}}^{(j)}(t)=\mathbf{{L}}^{(j)}w^{(j)}(t)$, is a vector-valued Gaussian process having spectral density $\mathbf{{Q}}_{c}^{(j)}=\mathbf{{L}}^{(j)}\sigma_{w}^{(j)}{(\mathbf{{L}}^{(j)})}^{T}$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$. Equation 45 can be written using block matrices as follows:

[TABLE]

where

[TABLE]

In shorthand notation, the augmented SSM as shown in Equation 46 can be represented by

[TABLE]

Here $\operatorname{\bm{{x}}}^{a}_{c}(t)\in\operatorname{\mathbb{R}}^{n_{a}}$ is the augmented state vector, $\tilde{\operatorname{\bm{{w}}}}(t)\in\operatorname{\mathbb{R}}^{n_{a}}$ is a concatenated vector-valued Gaussian process with spectral density $\mathbf{{Q}}_{c}$, as shown

[TABLE]

and matrices $\operatorname{\mathbf{{F}}}_{ac}\in\operatorname{\mathbb{R}}^{n_{a}\times n_{a}}$, $\mathbf{{L}}_{ac}\in\operatorname{\mathbb{R}}^{n_{a}\times\operatorname{\mathit{n}_{\mathit{f}}}}$ and $\mathbf{{H}}_{ac}\in\operatorname{\mathbb{R}}^{\operatorname{\mathit{n}_{\mathit{y}}}\times n_{a}}$ where $n_{a}=\operatorname{\mathit{n}_{\mathit{x}}}+m_{1}+\ldots+m_{\operatorname{\mathit{n}_{\mathit{f}}}}$.

4.3 Joint posterior inference of inputs and states

The discrete-time form of the augmented SSM (Equation 48) can be obtained as

[TABLE]

where the state-space matrices $\operatorname{\mathbf{{F}}}_{ad}=e^{\operatorname{\mathbf{{F}}}_{ac}\Delta t}$ and $\mathbf{{H}}_{ad}=\mathbf{{H}}_{ac}$. $\tilde{\operatorname{\bm{{w}}}}_{k-1}$ is a zero-mean Gaussian white noise vector representing the discrete-time form of $\tilde{\operatorname{\bm{{w}}}}(t)$ whose covariance is given by

[TABLE]

where $\mathbf{{\Psi}}_{a}(\tau)=e^{\operatorname{\mathbf{{F}}}_{ac}\tau}$ is the matrix exponential of the state-transition matrix. The integral in Equation 51 is solved using matrix fraction decomposition (see [43] for implementation details). Zero-mean Gaussian white noise $\operatorname{\bm{{w}}}_{k-1}$ with covariance $\operatorname{\mathbf{{Q}}}^{x}$ can be added to the process model in Equation 50 to account for unmodelled dynamics. The resulting modified form of Equation 50 can be written as

[TABLE]

where $\tilde{\operatorname{\bm{{w}}}}^{a}_{k-1}$ is the modified white noise vector with modified covariance $\operatorname{\mathbf{{Q}}}^{a}$ as defined below:

[TABLE]

$\operatorname{\bm{{w}}}_{k-1}$ and $\operatorname{\bm{{v}}}_{k}$ are assumed to be uncorrelated to each other and their joint covariance is expressed through Equation 11.

The posterior distribution of the states and the forces $p(\operatorname{\bm{{x}}}^{a}_{k}|\operatorname{\bm{{y}}}_{1:N},\operatorname{\bm{{\theta}}})=\mathcal{N}(\bm{{m}}^{a}_{k|N},\operatorname{\mathbf{{P}}}^{a}_{k|N})$ can be estimated with the classical Kalman filter (and smoother). Here, $\operatorname{\bm{{\theta}}}$ is a vector comprising all covariance function hyperparameters i.e. $\operatorname{\bm{{\theta}}}=\{\operatorname{\bm{{\theta}}}_{j}\}_{j=1}^{\operatorname{\mathit{n}_{\mathit{f}}}}$. The estimation is started from a Gaussian noise prior $\mathcal{N}(\bm{{m}}^{a}_{0|0},\operatorname{\mathbf{{P}}}^{a}_{0|0})$ where the covariance matrix $\operatorname{\mathbf{{P}}}^{a}_{0|0}$ has the following block diagonal form

[TABLE]

$\operatorname{\mathbf{{P}}}^{x}_{0|0}$ is the initial covariance matrix for the non-augmented state vector $\operatorname{\bm{{x}}}_{k}$ chosen according to some prior knowledge. The covariance matrix $\operatorname{\mathbf{{P}}}^{(j)}_{0|0}$ for the $j$th latent force is set equal to the steady state covariance matrix obtained using Equation 37, i.e. $\operatorname{\mathbf{{P}}}^{(j)}_{0|0}=\operatorname{\mathbf{{P}}}^{(j)}_{\infty}$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$.

4.4 Hyperparameter optimization for GP covariance functions

The structural model parameter matrices and that of the covariance functions corresponding to latent forces are transformed into parameters of the augmented GPLFM state-space model. As already stated in Section 2.2, the knowledge of the structural system parameters is assumed to be known a priori (either from FE model or from system identification), and therefore the augmented state estimation results will depend only on the parameters of the chosen covariance functions for the latent force components. In general, a parametric family of covariance functions is chosen and the hyperparameters are optimized based on the measurement data. Typical hyperparameters include lengthscale and signal variance for a standard family of covariance functions. The hyperparameters can be estimated in different ways, including maximization of marginal likelihood [21], maximum a posteriori [44], and Markov Chain Monte Carlo methods [45]. In this study, the optimized hyperparameters are obtained by maximizing the likelihood function based on the measurements. Maximum likelihood estimates of the hyperparameters (i.e. signal variance and lengthscale) of the covariance function(s) can be obtained by minimizing the negative log-likelihood (or maximizing the log-likelihood) of the measurements as follows:

[TABLE]

Using Kalman filter recursions (see B), the probabilities can be obtained as $p(\operatorname{\bm{{y}}}_{0}|\operatorname{\bm{{\theta}}})=\mathcal{N}\left(\operatorname{\bm{{y}}}_{0}|\mathbf{{H}}_{ad}(\operatorname{\bm{{\theta}}})\operatorname{\bm{{m}}}^{a}_{0|0}\right)$ and p(\operatorname{\bm{{y}}}_{k}|\operatorname{\bm{{y}}}_{0:k-1},\operatorname{\bm{{\theta}}})=\mathcal{N}\left(\operatorname{\bm{{y}}}_{k}\big{|}\mathbf{{H}}_{ad}(\operatorname{\bm{{\theta}}})\hat{\operatorname{\bm{{m}}}}^{a}_{k|k-1}\right) where $\operatorname{\bm{{m}}}^{a}_{0|0}$ is the estimate of the initial augmented state vector $\operatorname{\bm{{x}}}^{a}_{0}$ and $\hat{\operatorname{\bm{{m}}}}^{a}_{k|k-1}$ is the $k$th predicted augmented state vector. The minimization of the negative log-likelihood of the measurements can then be expressed in terms of the innovations $\bm{{e}}_{k}=\left(\operatorname{\bm{{y}}}_{k}-\mathbf{{H}}_{ad}(\operatorname{\bm{{\theta}}})\hat{\operatorname{\bm{{m}}}}^{a}_{k|k-1}\right)$ and the innovation covariance $\mathbf{{S}}_{k}$ (see [3] for more details) as:

[TABLE]

The expressions for $\bm{{e}}_{k}$ and $\mathbf{{S}}_{k}$ are provided in Equations 77 and 78 respectively. The minimization can be done using optimization tools such as MATLAB’s built-in functions fminunc or fmincon. It is noteworthy to mention that maximum likelihood optimization may get stuck in a local minimum; to avoid this one may need to start the optimization from different initial points. An algorithm depicting the steps involved in the proposed algorithm is shown in Algorithm 1.

4.5 Some properties of GPLFM

Property 1: The augmented state-space model used in Kalman filter-based approaches [15, 16, 19] can be considered a special case of the proposed GPLFM.

Proof.

The continuous-time SSM of the augmented Kalman filter formulation with noisy discrete measurements can be written as

[TABLE]

where $\tilde{\operatorname{\bm{{w}}}}^{f}(t)$ is a vector Gaussian process driving the evolution of the force $\operatorname{\bm{{f}}}(t)$. Comparing the above equation with the augmented SSM of GPLFM in Equation 46, one can deduce that the augmented SSM of GPLFM reduces to that of AKF in Equation 57 if $\operatorname{\bm{{z}}}(t)=\operatorname{\bm{{f}}}(t)$ or equivalently, $\operatorname{\mathbf{{F}}}^{*}=\operatorname{\mathbf{{0}}}_{\operatorname{\mathit{n}_{\mathit{f}}}\times\operatorname{\mathit{n}_{\mathit{f}}}}$ (due to $\operatorname{\mathbf{{F}}}^{(j)}=0$), $\operatorname{\mathbf{{B}}_{c}^{*}}=\operatorname{\mathbf{{B}}_{c}}$ and $\operatorname{\mathbf{{J}}_{c}^{*}}=\operatorname{\mathbf{{J}}_{c}}$ (due to $\mathbf{{H}}^{(j)}=1\;\forall\;j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$). Hence, the state-space model of AKF [15] can be considered as a special case of GPLFM. The state-space formulations in [16, 19] are similar to the augmented formulation in [15] and hence can be also considered special cases of GPLFM. ∎

Property 2: Assume that $(\operatorname{\mathbf{{A}}_{c}},\operatorname{\mathbf{{G}}_{c}})$ is observable, and that each latent force component $\operatorname{\bm{{z}}}^{(j)}(t),\;j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$ has an exponentially stable state-space representation. Then the augmented system $(\operatorname{\mathbf{{F}}}_{ac},\mathbf{{H}}_{ac})$ is detectable.

Proof.

Since the state-space representation $\left(\operatorname{\mathbf{{F}}}^{(j)},\mathbf{{H}}^{(j)}\right)$ of the $j$th latent force component $\operatorname{\bm{{z}}}^{(j)}(t)$ is exponentially stable, the combined state-space representation of all latent force components is also exponentially stable. The detectability of the full system $\left(\operatorname{\mathbf{{F}}}_{ac},\mathbf{{H}}_{ac}\right)$ can be determined using the Popov-Belevitch-Hautus (PBH) criterion ([46], which states that the full system is detectable if and only if the PBH matrix

[TABLE]

is a full column-rank for all $s\in\mathbb{C}$ or the undetectable modes are stable. This criteria has been used to show un-observability of the augmented state-space model given in Equation 57 when only acceleration measurements are used (see Equation (21) in [16]). In the case of GPLFM, expanding Equation 58 in terms of block matrices

[TABLE]

it can be seen that the state-space representation of GPLFM has full column-rank $\forall\;s\in\mathbb{C}$ due to $(\operatorname{\mathbf{{A}}_{c}},\operatorname{\mathbf{{G}}_{c}})$ being observable and $(\operatorname{\mathbf{{F}}}^{(j)},\mathbf{{H}}^{(j)})$ being stable for all $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$. Therefore, the augmented system of GPLFM is detectable with all types of measurements. ∎

Property 3: Assume that $\operatorname{\mathbf{{A}}_{c}}$ and $\operatorname{\mathbf{{F}}}^{(j)}$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$ are stable, then irrespective of the type of measurements (e.g. velocity, accelerations) used, the augmented system $(\operatorname{\mathbf{{F}}}_{ac},\mathbf{{H}}_{ac})$ has no marginally stable transmission zeros and admits stable inversion.

Proof.

The proof presented here uses the Moylan’s algorithm [47] for continuous-time linear systems. Consider the following matrix $\mathbf{{U}}(s)$

[TABLE]

The dimensions of $\operatorname{\mathbf{{A}}_{c}}$ are $\operatorname{\mathit{n}_{\mathit{x}}}\times\operatorname{\mathit{n}_{\mathit{x}}}$, $\operatorname{\mathbf{{B}}_{c}^{*}}$ is $\operatorname{\mathit{n}_{\mathit{x}}}\times M$, $\operatorname{\mathbf{{F}}}^{*}$ is $M\times M$, $\operatorname{\mathbf{{G}}_{c}}$ is $\operatorname{\mathit{n}_{\mathit{y}}}\times\operatorname{\mathit{n}_{\mathit{x}}}$ and $\operatorname{\mathbf{{J}}_{c}^{*}}$ is $\operatorname{\mathit{n}_{\mathit{y}}}\times M$. Note $M=m_{1}+\ldots+m_{\operatorname{\mathit{n}_{\mathit{f}}}}$, where $m_{j}$ is the size of $\operatorname{\mathbf{{F}}}^{(j)}$. The dimensions of $\mathbf{{U}}$ are $(\operatorname{\mathit{n}_{\mathit{x}}}+M+\operatorname{\mathit{n}_{\mathit{y}}})\times(\operatorname{\mathit{n}_{\mathit{x}}}+M)$, and therefore the maximum rank of $\mathbf{{U}}$ can be $(\operatorname{\mathit{n}_{\mathit{x}}}+M)$. Using the Schur complement of $(\operatorname{\mathbf{{A}}_{c}}-s\operatorname{\mathbf{{I}}})$, one can equivalently write

[TABLE]

Evaluating the rank of the above system at $s=0$, it can be seen that

[TABLE]

Since $\operatorname{\mathbf{{F}}}^{(j)}$, $j=1,\ldots,\operatorname{\mathit{n}_{\mathit{f}}}$ are assumed stable, $\operatorname{\mathbf{{F}}}^{*}$ is stable (implying none of the eigenvalues of $\operatorname{\mathbf{{F}}}^{*}$ are zero), and $\text{rank}(\mathbf{{F}}^{*})=M$. Therefore,

[TABLE]

irrespective of the rank of $\left(\operatorname{\mathbf{{J}}_{c}^{*}}-\operatorname{\mathbf{{G}}_{c}}\operatorname{\mathbf{{A}}_{c}}^{-1}\operatorname{\mathbf{{B}}_{c}^{*}}\right)$. Also, $\operatorname{\mathbf{{A}}_{c}}$ is assumed stable and has $\text{rank}(\operatorname{\mathbf{{A}}_{c}})=\operatorname{\mathit{n}_{\mathit{x}}}$, and hence

[TABLE]

This shows that the rank of the augmented GPLFM state-space model does not reduce at $s=0$. Therefore the GPLFM has no marginally stable transmission zeros. If the GP dynamics associated with the forces are removed, that is if $\operatorname{\mathbf{{F}}}^{(j)}=0$, then $\operatorname{\mathbf{{F}}}^{*}$ degenerates to $\operatorname{\mathbf{{0}}}_{\operatorname{\mathit{n}_{\mathit{f}}}\times\operatorname{\mathit{n}_{\mathit{f}}}}$, and $\operatorname{\mathbf{{B}}_{c}^{*}}$ and $\operatorname{\mathbf{{J}}_{c}^{*}}$ reduces to $\operatorname{\mathbf{{B}}_{c}}$ and $\operatorname{\mathbf{{J}}_{c}}$ respectively, and one reverts to the state-space model of AKF (see Property 1). The state-space model underlying the AKF has been shown to possess marginally stable transmission zeros when used with only acceleration/velocity measurements [48], which typically leads to drifts in force estimation. Since the GPLFM has no marginally stable transmission zeros, its inversion is stable and is robust to drift in force estimation. ∎

5 On the choice of covariance functions

In practice, for modelling the forces one has to first choose the covariance functions and then optimize the hyperparameters such that the mismatch between the predicted responses and the observed responses is minimized. In general, this choice is determined by the expected properties of the underlying forces. In this study, a valid choice of covariance function must satisfy these two criteria: (i) they must belong to the class of stationary covariance functions, and (ii) they must admit a stable state-space representation following the spectral factorization approach as described in Section 4.1. Several covariance functions exist that satisfy the two criteria, such as the Matérn class, the squared exponential, the exponential, the periodic and the rational quadratic covariance functions. An illustration of different covariance functions and their output processes are shown in Figure 2. An added advantage of GPLFM is that different individual covariance functions can be combined to create a rich set of covariance functions, that is, any sum and/or product of the above covariance functions gives a valid covariance function. For example, the sum of squared exponential and periodic covariance functions can be used to model forces that have both random and periodic components. Thus, the set of available covariance functions not only includes the list of individual covariance functions but also all valid combinations constructed by taking sum and/or product between them. However, this comes at the cost of increased complexity due to an increase in the number of tunable hyperparameters.

Typically, one would choose a covariance function (or a valid combination of covariance functions) that best fits some training data (i.e., measurements of the true force). In the absence of training data, a covariance function can be chosen based on expert guess or prior knowledge. In the absence of both, one would resort to popular choices such as the Matérn family of stationary covariance functions. The Matérn family of stationary covariance functions offers great flexibility in modelling different types of random processes while having only a few tunable hyperparameters. The exponential and the squared exponential covariance functions are also special cases of the Matérn class of covariance functions. Illustrations of different Matérn covariance functions are provided in Figure 3.

In the field of structural dynamics, the operational forces commonly sustained by structures can be modelled as random processes. For example, wind and seismic loads are typically non-stationary and narrow-banded random processes. Such non-stationary forces can be effectively modelled as GPs generated from stationary covariance functions. For example, rapidly varying (smoothly varying) random forces can be satisfactorily modelled by the GP with exponential (squared exponential) covariance functions. If in addition to randomness there are periodic trends, then a customized covariance function can be constructed by taking the sum/product of periodic and squared exponential covariance functions. However, care needs to be taken when modelling non-random loads that are characterized by abrupt changes or ‘kinks’ in their time history (e.g. impact loads, step loads). As stationary covariance functions are usually continuous and differentiable everywhere, capturing the kinks or abrupt changes using GPs with stationary covariance functions can be difficult. In such cases, one may be obliged to adopt the exponential covariance function (which is continuous everywhere but not differentiable) or a combination of the exponential covariance function with the squared exponential function as a useful choice.

6 Numerical studies

For investigation of the performance of the proposed GPLFM, a 10-storey building is chosen. Each floor of the building is represented by a dof leading to a 10-dof shear building model. The mass of each floor is assumed to be $200$\mathrm{k}\mathrm{g} and the inter-storey floor stiffness as $5\times 10^{5}$\mathrm{N}\mathrm{/}\mathrm{m}. The damping matrix matrix $\mathbf{{C}}$ is assumed to be Rayleigh distributed as $\mathbf{{C}}=0.1\mathbf{{M}}+0.0005\mathbf{{K}}$. The modal frequencies and damping ratios are provided in Table 1.

Only acceleration measurements are assumed available at the floor levels of the structure. For the numerical simulations in this section, the acceleration responses of the structure have been sampled at 100 $\mathrm{H}\mathrm{z}$ and contaminated with 10% white Gaussian noise to model measurement errors.

The performance of the proposed GPLFM method has been compared to a suite of Kalman filter-based joint input-state estimation methods, namely, the augmented Kalman filter (AKF) [15], the dual Kalman filter (DKF) [18], and the augmented Kalman filter with dummy measurements (AKFdm) [19]. All the Kalman filter based methods need the a priori knowledge of the mean and covariance of the initial state $\left(\operatorname{\bm{{m}}}^{x}_{0|0},\operatorname{\mathbf{{P}}}^{x}_{0|0}\right)$, the mean and covariance of the input $\left(\operatorname{\bm{{m}}}^{f}_{0|0},\operatorname{\mathbf{{P}}}^{f}_{0|0}\right)$, the state and input process noise covariance matrices $\left(\operatorname{\mathbf{{Q}}}^{x},\operatorname{\mathbf{{Q}}}^{f}\right)$, and the measurement noise covariance matrix $\operatorname{\mathbf{{R}}}$. In addition, AKFdm — an extension of AKF, which was proposed to alleviate low frequency drift in force estimates encountered in AKF when using only accelerations — needs the specification of a dummy displacement measurement covariance matrix $\operatorname{\mathbf{{R}}}_{dm}$. Mean values of the initial state $\operatorname{\bm{{m}}}^{x}_{0}$ and the input $\operatorname{\bm{{m}}}^{f}_{0}$ are typically assigned zero values for most purposes. The initial state covariance $\operatorname{\mathbf{{P}}}^{x}_{0}$ and the initial input covariance $\operatorname{\mathbf{{P}}}^{f}_{0}$ represent the uncertainty in the initial state and input respectively. A lower (or higher) covariance implies stronger (or weaker) belief in the assumed initial values. Process noise covariance matrix $\operatorname{\mathbf{{Q}}}^{x}$ signifies the uncertainty in the assumed dynamical model of the structure arising primarily due to modelling error, and lower (or higher) values imply higher (or lower) confidence in the dynamic model. In comparison to $\operatorname{\mathbf{{Q}}}^{x}$, the action of input process noise covariance matrix $\operatorname{\mathbf{{Q}}}^{f}$ is a bit different, in that it represents the covariance of a fictitious noise process through which the input estimate evolves and plays the role of a regularizer that strongly influences the performance of input and state estimation in AKF, DKF and AKFdm algorithms. A higher $\operatorname{\mathbf{{Q}}}^{f}$ will allow the input to assume larger values while varying rapidly, whereas a lower $\operatorname{\mathbf{{Q}}}^{f}$ will restrict the input to take relatively smaller values while varying more gradually. In this study, calibration of the values of $\operatorname{\mathbf{{Q}}}^{f}$ for use in AKF, DKF and AKFdm has been carried out mostly with the help of L-curve method along with some manual tuning. The L-curve is obtained by plotting the summed mean squared values of the innovation sequence i.e. $\sum\left\lVert\bm{{e}}_{k}\right\rVert$ against the values of $\left\lVert\operatorname{\mathbf{{Q}}}^{f}\right\rVert$, and the values of $\operatorname{\mathbf{{Q}}}^{f}$ corresponding to the “corner” of the L-curve are selected. The measurement noise covariance matrix $\operatorname{\mathbf{{R}}}$ represents the degree of uncertainty in the measurements due to measurement instrument errors; an estimate of $\operatorname{\mathbf{{R}}}$ can be obtained from the precision of the instrument prescribed by the manufacturer. The values of the dummy displacement measurement covariance matrix $\operatorname{\mathbf{{R}}}_{dm}$, needed for AKFdm, has been recommended [16] to be chosen with an order of magnitude higher than the covariance of displacement responses for stable estimation.

In GPLFM, one needs to provide $\operatorname{\bm{{m}}}^{x}_{0|0}$, $\operatorname{\mathbf{{P}}}^{x}_{0|0}$, $\operatorname{\mathbf{{Q}}}^{x}$, $\operatorname{\mathbf{{R}}}$, $\operatorname{\bm{{m}}}^{f}_{0|0}$ and the choice of the family of covariance function for latent forces prior to implementation. Unlike DKF, AKF and AKFdm where the tuning of $\operatorname{\mathbf{{Q}}}^{f}$ controls the accuracy of estimation results, the proposed GPLFM method does not need the specification of $\operatorname{\mathbf{{Q}}}^{f}$; instead in GPLFM the hyperparameters of chosen covariance functions dictate the performance of input and state estimation. The hyperparameters of the covariance function can be set based on expert knowledge or computed through maximum likelihood optimization based on measured data as described in Section 4.4. The latter has been adopted in the study to avoid the dependence on any prior guess or expert knowledge. Once the optimal hyperparameters are computed, the input and state estimation is performed using Kalman filter.

In this study, exponential covariance functions —a class of the popular Whittle-Matérn family of stationary covariance functions (see A) —have been used for GP modelling of the different excitations, unless explicitly mentioned otherwise. Exponential covariance functions are non-differentiable and the GPs generated using them are continuous everywhere but not differentiable, and are thus useful in modelling random (less smooth) phenomena such as Brownian motion. Since most excitations sustained by operational structures are random in nature, GPs with exponential covariance functions are employed for modelling the input excitations in this study.

For all numerical simulations in this section, the following values are assumed: $\bm{{m}}_{0|0}^{x}=\operatorname{\mathbf{{0}}}$, $\bm{{m}}_{0|0}^{f}=\operatorname{\mathbf{{0}}}$, $\operatorname{\mathbf{{P}}}^{x}_{0|0}=10^{-10}\times\operatorname{\mathbf{{I}}}$, $\operatorname{\mathbf{{Q}}}^{x}=10^{-10}\times\operatorname{\mathbf{{I}}}$ and $\operatorname{\mathbf{{R}}}=0.1\times\operatorname{\mathbf{{I}}}\;($\mathrm{m}\mathrm{/}\mathrm{s}^{2}$)^{2}$ for GPLFM, DKF, AKF and AKFdm. In what follows, the performance of joint input and state estimation is assessed using different types of excitations i.e. impact, harmonic, seismic and random.

6.1 Impact excitation

An impact excitation of magnitude $10^{4}$\mathrm{N} is applied at the topmost floor of the structure. The impact load starts at $t=3$\mathrm{s} of the simulation, ramps up linearly from 0 to $10^{4}$ in 0.05s, peaks at $t=3.05$\mathrm{s} and then ramps downs linearly from $10^{4}$ to 0 in another $0.05$\mathrm{s}. The following values of covariances are used in input and state estimation

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{4}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for DKF, AKF and AKFdm obtained using L-curve criterion as shown in Figure 4,

• •

$\operatorname{\mathbf{{R}}}_{dm}=0.05\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm obtained by setting the covariance roughly 10 times the square of the maximum absolute displacement value.

Two cases of measurements are considered:

1. (a)

All acceleration measurements

In this case, the state and input estimation is performed using acceleration measurements from all floor levels. The estimated acceleration, velocity and displacement time histories at the 5th floor (chosen as a representative floor) is shown in Figure 5. It is found that all the algorithms are able to accurately estimate the acceleration and the velocity states, however the displacement states estimated by AKF and DKF are affected by low frequency drifts. The low frequency drift that manifests in the AKF and the DKF estimates are due to accumulation of integration errors. The AKFdm which uses dummy displacement measurements is able to arrest the drift and provide stable displacement estimates. Both AKFdm and GPLFM provide reasonably accurate estimates of the displacements, with GPLFM performing the best among all the algorithms. Similar estimation results are obtained at all other floor levels. It can be seen that state estimates provided by AKFdm compares closely with that of AKF except the displacement estimates where drift errors due to AKF leads to large estimation error in estimates obtained by AKF.

The estimated input force histories are shown in Figure 6. It is seen that input force estimates using DKF and AKF are affected by low frequency drift, alike the displacement estimates. AKFdm force estimate shows some post peak oscillations but does well in restricting the drift in estimates. GPLFM seem to provide good force estimates with no drift and negligible oscillations. It can be seen that none of the estimates from any of these algorithms could achieve the exact peak of the actual impulsive force, with the GPLFM estimate furnishing the closest estimate to the true force.

2. (b)

*Single acceleration measurement

*Only a single acceleration at the 10th floor level is measured in this case. This corresponds to a collocated measurement case where the force and the measurement locations coincide. Collocated acceleration measurements are desirable in force estimation since the input force contributes directly to collocated acceleration measurements through the direct feedthrough term. The estimates of force time history are provided in Figure 7. The force estimates obtained using AKF and DKF show low frequency drifts as was encountered previously. However, for AKF, the estimate of the force at the peak improved in this case. The peak force estimated by AKF, DKF and AKFdm match closely to each other, although underestimated by around 40%. The peak force estimated by GPLFM appear slightly better than the estimated by DKF, AKF and AKFdm, however, the peak is underestimated by 30% when using a single acceleration measurement. A post peak trough is observed close to 3.5$\mathrm{s}$ in the force estimates from all algorithms. In the authors’ opinion, the appearance of the trough may have lead to a drop in the peak force estimated by these algorithms.

The acceleration, velocity and displacement state estimates at the 5th floor is shown in Figure 8. It is found that even with the use of a single collocated acceleration, the acceleration and velocity state estimates obtained by all the algorithms at other floors still show good accuracy. For the displacement state, however, the GPLFM provides a better estimate compared to the estimates from other algorithms.

It should be mentioned here that the attempt to estimate an impact loading by a Gaussian process with exponential covariance function is a scenario of mis-specified covariance function model. The mis-specification arises because the impact excitation is being modelled here by a Gaussian process with a continuous exponential covariance function, although an impact is best represented by a discontinuous Dirac-delta function, atleast in theory. Such mis-specification should be avoided if possible, however, for practical reasons such as convenience of using standard forms of covariance function or vague prior information, one inevitably ends up in a situation which resembles some level of mis-specification [21].

6.2 Harmonic excitation

In this section, the performances of GPLFM, DKF, AKF and AKFdm are assessed using harmonic excitation. A constant amplitude sinusoidal excitation of the form $f(t)=A_{h}\sin(2\pi\lambda t)$ is applied to the topmost floor of the shear building, where $A_{h}=100$\mathrm{N}, $\lambda=1$\mathrm{H}\mathrm{z}. The following values are used prior to running the algorithms:

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{4}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for DKF, AKF and AKFdm set using L-curve criterion,

• •

$\operatorname{\mathbf{{R}}}_{dm}=10^{-2}\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at 100 times the absolute maximum of true displacement history,

• •

Matérn covariance function with $p=2$ (refer to Equation 73) is chosen for GPLFM

For the sinusoidal case, a Matérn covariance function with $p=2$ is used in place of exponential covariance function as it was found difficult to obtain good force estimates with a exponential covariance function. Ideally, a periodic covariance function is most appropriate for estimating a periodic force history, however, a more convenient Matérn covariance function with $p=2$ having the property of being twice differentiable is used in this case. Once again this partly resembles a mis-specified covariance function.

Using a single collocated acceleration measurement at the topmost floor, the force history is estimated as shown in Figure 9. It is found that GPLFM is able to provide very accurate estimate of the sinusoidal force history. Furthermore, it can be seen that DKF performed better than AKF in estimating the magnitude of the force, however, both their estimates are affected by drift. AKFdm force estimate tracks the true force in an unbiased manner, although some low frequency oscillations ride the force estimate. Reducing the value of $\operatorname{\mathbf{{R}}}_{dm}$ removed the low frequency oscillations in the force estimate from AKFdm, however, that occured at the cost of poor estimation of the force magnitude. The state estimates obtained by the algorithms at the representative 5th floor compare reasonably well with respect to each other as seen in Figure 10. The drift in the estimates of DKF and AKF in this case is found to be comparably low. The estimation case with all accelerations produced similar results and is thus not included here.

6.3 Seismic excitation

To assess the performance of different algorithms under seismic excitation, the El Centro earthquake record (Imperial Valley Station, May 18, 1940) is used as ground acceleration input to the 10-storey structure. Only the case of a single acceleration measurement at the 10th floor is considered. It was found that using multiple floor accelerations gave similar results and hence not discussed. It is important to note that accelerometers measure only absolute accelerations, and in this case, when no other external force is acting on the structure, measuring absolute accelerations reduces the direct feedthrough term in the observation equation to zero (refer to derivation in D). When the direct feedthrough term becomes zero, the DKF algorithm fails to estimate the input since the feedthrough matrix $\mathbf{{J}}$ in Equation 10 becomes zero and this leads to zero Kalman gain for the input, thereby degenerating the input update stage in DKF implementation. Thus, for this case, a comparison will be drawn among estimates from GPLFM, AKF and AKFdm. The following values are adopted prior to running the algorithms:

• •

$\operatorname{\mathbf{{Q}}}^{f}=0.1\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for AKF and AKFdm obtained using L-curve criterion,

• •

$\operatorname{\mathbf{{R}}}_{dm}=\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at roughly 10 times the square of absolute maximum of true displacement history,

The estimated ground acceleration histories from GPLFM, AKF and AKFdm are shown in Figure 11. The GPLFM estimate is able to track the ground acceleration time history quite accurately. Both AKF and AKFdm also demonstrate good tracking of the ground acceleration although with a bit of noise in the estimates.

The estimated acceleration, velocity and displacement states at the 5th floor are shown in Figure 12. The state estimates from all the three algorithms at the 5th floor seem reasonably accurate with GPLFM furnishing the closest estimate to the true states.

6.4 Random excitation

Most operational structures are subjected to random forces (such as wind, traffic, etc.) that can be usefully modelled by random white noise (or filtered white noise) excitation. In this subsection, the performance of GPLFM, DKF, AKF and AKFdm in input and state estimation is assessed under the application of random excitation to the 10-storey shear structure. Two scenarios of random excitation are considered: (I) a single random excitation applied at the 10th floor of the structure, and (II) multiple random excitation applied to all floors of the structure.

1. (I)

*Single random excitation

*A zero mean Gaussian white noise excitation with standard deviation $1$\mathrm{k}\mathrm{N}$$ is applied to the 10th floor. The following values of covariances are used:

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{4}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for DKF, AKF and AKFdm obtained using L-curve criterion and some manual tuning,

• •

$\operatorname{\mathbf{{R}}}_{dm}=0.05\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at 100 times the square of absolute maximum of true displacement history

First consider the case where accelerations at all floor levels are measured. Using all acceleration measurements, the estimated force histories from GPLFM, DKF, AKF and AKFdm are shown in Figure 13 and the estimated states for the 5th floor are shown in Figure 14. It is found that GPLFM estimate is able to provide a very good tracking of the true force history compared to the estimate from other algorithms. As regards the state estimates, both GPLFM and AKFdm provide good estimates of all states while AKF and DKF estimates for displacement states suffer from drift.

Next, the case of a single collocated acceleration measured at the 10th floor is considered. The force and state estimation results from GPLFM, DKF, AKF and AKFdm obtained for this measurement case are shown in Figure 15 and 16 respectively. It is found that estimation results from using a single collocated acceleration are still reasonable but not as accurate when compared to the results from using all accelerations.

Now, a case where accelerations are measured at floor levels 1, 3, 5, 7 and 9 is considered. This represents a scenario of non-collocated measurements, where the locations of measurements do not coincide with the location(s) of the applied force(s). The force and state estimation results for this case are shown in Figure 17 and 18 respectively. It is found that input force estimation using non-collocated observations proves to be a challenging case. The DKF fails to estimate the forces due to the direct feedthrough matrix being zero. The force estimates from GPLFM, AKF and AKFdm seem to somewhat follow the true force history but still are far from accurate, however, among them the GPLFM estimate was found to have the least mean square error. On the other hand, the state estimation results from all the algorithms are found to be reasonably accurate except the state estimates from DKF.

2. (II)

*Multiple random excitation

*Here, a scenario of multiple random forces exciting the structure at all floor levels is considered. Random excitations having zero mean and standard deviation $0.1$\mathrm{k}\mathrm{N}$$ are applied at all floor levels, and acceleration measurements are assumed to be available at floors 1, 3, 5, 7 and 9. This represents a case where (a) the number of forces to be estimated are greater than the number of observations (10 forces against 5 observations), and (b) the observations at floor levels 1, 3, 5, 7 and 9 serve as both collocated observations and non-collocated observations; collocated for forces 1, 3, 5, 7, 9 and non-collocated for forces 2, 4, 6, 8, 10.

For this the following values of covariances are chosen:

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{3}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for DKF and AKF chosen using L-curve criteron and some manual tuning,

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{4}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ AKFdm chosen after some manual tuning,

• •

$\operatorname{\mathbf{{R}}}_{dm}=10^{-2}\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at approximately two times the square of the absolute maximum of true displacement.

The estimated force time histories for two representative floors, 8th floor (non-collocated measurement) and 9th floor (collocated measurement), are illustrated in Figure 19. It can be seen that for the 8th floor, none of the algorithms are able to estimate the true force history correctly. Both the GPLFM and the AKF seem to provide no tracking of the true force at the 8th floor, besides AKF suffering from drift. The DKF is unable to estimate the forces at non-collocated measurement locations due to direct transmission term becoming zero leading to degeneracy in the input update stage of DKF. The AKFdm estimate of the force at the 8th floor seem to produce larger magnitudes but no accuracy in tracking the true force. Comparing this with the force estimates (in Figure 19) obtained at the 9th floor, it can be seen that GPLFM estimates of force history is quite accurate while AKF and DKF estimates for the collocated forces are able to somewhat trace the pattern of the force history albeit with significant drift.

The displacement and velocity state estimates at the 8th and 9th floors are shown in Figure 20. Both DKF and AKF are able to provide accurate velocity estimates alike GPLFM but suffer from drift in displacement estimates although they seem to follow the pattern of the true displacement state quite accurately. The GPLFM estimates for displacements and velocities are quite accurate both for the case of collocated and non-collocated measurements.

7 Application: A 76-storey ASCE benchmark building

An ASCE benchmark 76-storey office tower [49] is considered here for application. The 76-storey office tower was proposed for the city of Melbourne, Australia. The building plan and elevation view have been provided in Figure 21.

The 76-story tall building is modelled as a vertical cantilever beam. A finite element model is constructed by considering the portion of the building between two adjacent floors as a classical beam element of uniform thickness, leading to 76 translational and 76 rotational degrees of freedom. Then, all the 76 rotational degrees of freedom have been removed by the static condensation. This resulted in 76 degrees of freedom, representing the translation of each floor in the lateral direction. The first five natural frequencies are 0.16, 0.765, 1.992, 3.790 and 6.395 $\mathrm{H}\mathrm{z}$. Damping ratios for the first five modes are assumed to be 1%. This model, having $76\times 76$ mass, damping and stiffness matrices, is referred to as the 76-dof ‘true’ model. For simplification of numerical computation, the 76-dof true model was reduced to a 23-dof system such that the first 46 ($=2\times 23$) complex modes of the 76-DOF model are retained. The 23-dof system correspond to the translation at floors 3, 6, 10, 13, 16, 20, 23, 26, 30, 33, 36, 40, 43, 46, 50, 53, 56, 60, 63, 66, 70, 73 and 76, and is referred to as the 23-dof reduced-order model. Responses at nine floor levels are assumed to be available — floor levels 1, 30, 50, 55, 60, 65, 70, 75 and 76. All relevant details and MATLAB files for the structure are provided in the website [50] for ASCE benchmark buildings.

Two separate scenarios of excitation, seismic and wind, are considered in this application. For the case of seismic excitation, the 76-dof true model is subjected to El Centro earthquake ground motion, and absolute acceleration responses, sampled at 100$\mathrm{H}\mathrm{z}$, are generated at floor levels 1, 50 and 70, for a duration of 50$\mathrm{s}$. The generated acceleration responses are then contaminated with 10% measurement noise. These noisy measurements along with the 23-dof reduced-order model are used by GPLFM, AKF and AKFdm for input and state estimation. The following values are used for this scenario:

• •

$\operatorname{\mathbf{{P}}}^{x}_{0|0}=10^{-10}\times\operatorname{\mathbf{{I}}}$ and $\operatorname{\mathbf{{Q}}}^{x}=10^{-10}\times\operatorname{\mathbf{{I}}}$ for GPLFM, AKF and AKFdm,

• •

$\operatorname{\mathbf{{R}}}=0.5\times\operatorname{\mathbf{{R}}}_{true}$ for GPLFM, AKF and AKFdm, where $\operatorname{\mathbf{{R}}}_{true}$ is the true measurement noise covariance with the $i$th diagonal element obtained as $\operatorname{\mathbf{{R}}}_{true}(i,i)=0.01*\text{Var}(\operatorname{\bm{{y}}}_{i})$

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{3}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for AKF and AKFdm obtained using L-curve criterion and some manual tuning,

• •

$\operatorname{\mathbf{{R}}}_{dm}=10\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at approximately 10 times the square of absolute maximum of true displacement history,

• •

An exponential covariance function (refer to Equation 71) is used for GPLFM

The estimated ground input time histories are shown in Figure 22. The estimated acceleration, velocity and displacement states at floor level 30, chosen as the representative floor, are shown in Figure 23. It is found that all the three algorithms i.e. GPLFM, AKF and AKFdm are able to track the ground motion input with good accuracy, with estimates from AKF and AKFdm showing slight overshoots. The state estimates obtained by GPLFM and AKF also match the true states quite accurately while the displacement estimates from AKFdm are under-estimated.

Next, the 76-storey model is subjected to wind excitation at the topmost floor i.e. 76th floor. For wind excitation, 200$\mathrm{s}$ of wind loading history was obtained from the 76th column of wind_cross.mat data file and then multiplied with a factor of 100; the multiplying factor was chosen such that the displacement obtained at the 76th floor was less than 0.5$\mathrm{m}$. Two cases of measurements are considered: (a) three accelerations measured at floor levels 1, 50 and 70, and (b) one displacement measured at 1st floor in addition to three accelerations measured at floor levels 1, 50 and 70. Note that these measurements now resemble a non-collocated measurement scenario. The GPLFM, AKF and AKFdm are used for input and state estimation for these cases; results from DKF are omitted as it DKF has been shown to suffer from degeneracy when non-collocated acceleration measurements are used. The following modified values are used in estimation:

• •

$\operatorname{\mathbf{{P}}}^{x}_{0}=10^{-10}\times\operatorname{\mathbf{{I}}}$ and $\operatorname{\mathbf{{Q}}}^{x}=10^{-10}\times\operatorname{\mathbf{{I}}}$ for GPLFM, AKF and AKFdm,

• •

$\operatorname{\mathbf{{R}}}=0.5\times\operatorname{\mathbf{{R}}}_{true}$ for GPLFM, AKF and AKFdm, where $\operatorname{\mathbf{{R}}}_{true}$ is the true measurement noise covariance with the $i$th diagonal element obtained as $\operatorname{\mathbf{{R}}}_{true}(i,i)=0.01*\text{Var}(\operatorname{\bm{{y}}}_{i})$

• •

$\operatorname{\mathbf{{Q}}}^{f}=10^{10}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for AKF and $\operatorname{\mathbf{{Q}}}^{f}=2\times 10^{7}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ for AKFdm obtained with manual tuning,

• •

$\operatorname{\mathbf{{R}}}_{dm}=3\times\operatorname{\mathbf{{I}}}\;$\mathrm{m}$^{2}$ for AKFdm set at approximately 10 times the square of absolute maximum of true displacement history.

• •

An exponential covariance function is used for GPLFM

1. (a)

*Three acceleration measurements

*The estimated displacement, velocity and acceleration states at the 30th floor level are shown in Figure 24. It is seen that the state estimates from GPLFM and AKF are reasonably close to the true value, however, the velocity and displacement estimates from AKFdm do not match the true states that closely.

From the estimated force time histories shown in Figure 25, it is apparent that none of the algorithms are able to provide proper tracking of the true input force, nevertheless, the state estimates obtained from GPLFM and AKF are still reasonable. Using nine more acceleration measurements did not help much in improving the force estimates. To improve the force estimation, the use of an additional displacement measurement is considered next. 2. (b)

*One displacement and three acceleration measurements

*A displacement measurement at the 1st floor is added to three accelerations measured at floors 1, 50 and 70 to improve upon the input force estimate from the previous case. A readjusted value of $\operatorname{\mathbf{{Q}}}^{f}=10^{2}\times\operatorname{\mathbf{{I}}}\;$\mathrm{N}$^{2}$ is now used for AKF and AKFdm. The estimated forces time histories are shown in Figure 26(a). The force estimates obtained from GPLFM, AKF and AKFdm are very close to each other and seem to track the force history with a constant delay of around 2$\mathrm{s}$. This delayed behavior can sometimes be expected in force estimation approach as reported in [51, 52]. The delayed behavior can be removed by applying a Kalman smoother as shown in Figure 26(b).

The acceleration, velocity and displacement state estimates for the 30th floor show good agreement with the true values as depicted in Figure 27.

8 Conclusions

In this paper, a novel implementation of GPLFM for joint input-state estimation of linear time-invariant structural systems with unknown inputs is proposed. The unknown inputs to a linear structural system are modelled as GPs with specified covariance functions. The GPLFM is then formulated into an augmented state-space model with additional states representing the latent force components. The GPLFM state-space model so obtained is shown to be a generalization of the augmented state-space models used in popular Kalman filter-based joint input-state estimation algorithms. Moreover, the GPLFM formulation is analytically proven to be observable and robust to drift in force estimation even when used with only acceleration measurements. The behavior of the Gaussian process latent force components are determined by the hyperparameters of the chosen covariance functions. In this study, these hyperparameters are tuned using maximum likelihood optimization based on the observed data avoiding the need to set them based on expert guess.

The effectiveness and performance of GPLFM is demonstrated through several numerical studies using different excitations and measurement scenarios on a 10-dof shear building as well as a numerical model of the 76-dof ASCE benchmark office tower. The performance of GPLFM is compared to that of DKF, AKF and AKFdm. It is found from the numerical studies that GPLFM is able to deliver very accurate estimates of inputs and states in most cases. Moreover, GPLFM does not suffer from drift in force and/or displacement estimates, unlike AKF and DKF. The choice of covariance function, however, plays an important role in the accuracy of the input estimation and relies on the user’s ability to encode prior information of the behavior of the latent forces such as discontinuity, periodicity, etc. Finally, it should be mentioned that similar to Kalman filter-based joint input-state estimation approaches such as AKF, DKF and AKFdm, GPLFM is amenable to online implementation once the hyperparameters of the covariance function are properly calibrated.

9 Acknowledgments

Funding for this study through the Natural Sciences Engineering Research Council of Canada (Canada Research Chairs program) is gratefully acknowledged by the authors.

Appendix A Matérn class of covariance functions

A commonly used class of covariance functions is the Matérn family [21, 53, 54]

[TABLE]

where $\alpha^{2}$ and $l$ are positive hyperparameters that denote the signal variance and lengthscale respectively, and $\nu$ is a parameter controlling the smoothness of the Gaussian process. In this paper, the set of hyperparameters of a covariance function is denoted by $\operatorname{\bm{{\theta}}}=\{\alpha^{2},l\}$. These hyperparameters control the variability of the underlying Gaussian process. $\Gamma(\nu)$ is the Gamma function, $K_{\nu}$ is a modified Bessel function of second kind [55]. For one-dimensional processes, the spectral density of the Matérn covariance function (Equation 65) is

[TABLE]

where $\lambda=\sqrt{2\nu}/l$, and in this work $\nu=p+1/2$, where $p$ is a non-negative integer. Thus,

[TABLE]

This has a proper rational fraction form and can be spectrally factorized as

[TABLE]

from which a stable transfer function can be obtained as

[TABLE]

The corresponding spectral density of the Gaussian white noise process $w(t)$ is

[TABLE]

1. (a)

Setting $p=0$, the exponential covariance function is obtained with $\lambda=1/l$, $\sigma_{w}=2\alpha^{2}/l$ and the corresponding LTI model matrices in Equation 34 read as:

[TABLE] 2. (b)

Setting $p=1$, the Matérn class with $\nu=3/2$ is obtained with $\lambda=\sqrt{3}/l$, $\sigma_{w}=12\sqrt{3}\alpha^{2}/l^{3}$ and the corresponding LTI model matrices in Equation 34 read as:

[TABLE] 3. (c)

Setting $p=2$, the Matérn class with $\nu=5/2$ is obtained with $\lambda=\sqrt{5}/l$, $\sigma_{w}=400\sqrt{5}\alpha^{2}/3l^{5}$ and the corresponding LTI model matrices in Equation 34 read as:

[TABLE]

Appendix B Kalman filter and smoother

Consider a discrete-time time-invariant linear system with additive Gaussian noise

[TABLE]

where $\operatorname{\bm{{x}}}_{0}\sim\mathcal{N}(\hat{\operatorname{\bm{{m}}}}_{0|0},\operatorname{\mathbf{{P}}}_{0|0})$ is the initial state distribution, $\operatorname{\bm{{w}}}_{k-1}\sim\mathcal{N}(\operatorname{\mathbf{{0}}},\operatorname{\mathbf{{Q}}})$ and $\operatorname{\bm{{v}}}_{k}\sim\mathcal{N}(\operatorname{\mathbf{{0}}},\operatorname{\mathbf{{R}}})$ are the process noise and the measurement noise distribution respectively, with zero cross-correlation between $\operatorname{\bm{{w}}}_{k-1}$ and $\operatorname{\bm{{v}}}_{k}$. For such a system, the Kalman filter [56] computes the mean and the covariance of posterior filtering distribution sequentially using the following steps:

[TABLE]

for $k=1,\ldots,N$. Here $\hat{\operatorname{\bm{{m}}}}_{k|k-1}$ and $\hat{\operatorname{\bm{{m}}}}_{k|k}$ represent the $k$th predicted and filtered state estimate, respectively, and, $\operatorname{\mathbf{{P}}}_{k|k-1}$ and $\operatorname{\mathbf{{P}}}_{k|k}$ denote the $k$th predicted and filtered state error covariance matrices respectively. The recursion in Kalman filter is started from initial state $\hat{\operatorname{\bm{{m}}}}_{0|0}$ and covariance $\operatorname{\mathbf{{P}}}_{0|0}$. Following filtering step, the (fixed interval) smoothing step is given by the Rausch-Tung-Striebel (RTS) smoother [57]:

[TABLE]

for $k=N-1,\ldots,1$.

Appendix C Derivation of $\mathbf{{K}}_{xf}(t,t^{\prime})$

A complete derivation of the expression in Equation 30c is provided here.

[TABLE]

where in the second last line, the initial state $\operatorname{\bm{{x}}}_{0}$ is assumed to be independent of the force applied $\operatorname{\bm{{f}}}(t)$ leading to $\mathbb{E}\left[\operatorname{\bm{{x}}}_{0}\operatorname{\bm{{f}}}(t^{\prime})^{T}\right]=\mathbb{E}\left[\operatorname{\bm{{x}}}_{0}\right]\mathbb{E}\left[\operatorname{\bm{{f}}}(t^{\prime})^{T}\right]$ and then using $\mathbb{E}\left[\operatorname{\bm{{f}}}(t)\right]=\operatorname{\mathbf{{0}}}$ (from the definition of GP over $\operatorname{\bm{{f}}}(t)$ from Eqn 22), we get $\mathbb{E}\left[\operatorname{\bm{{x}}}_{0}\right]\mathbb{E}\left[\operatorname{\bm{{f}}}(t^{\prime})^{T}\right]=\operatorname{\mathbf{{0}}}$.

Appendix D Direct-feedthrough matrix in case of seismic excitation

When a structural system is subjected to earthquake ground motion, there is no external physical force acting on the structure. The equation of motion of the structure, when considering absolute motion $\bm{{u}}_{a}(t)$ of the structure, is given by

[TABLE]

where $\bm{{u}}_{g}(t)\in\operatorname{\mathbb{R}}^{n_{g}}$ is the time-history of ground motion caused due to earthquake. Often for convenience, one rewrites Equation 85 in terms of relative motion, $\bm{{u}}(t)=\bm{{u}}_{a}(t)-\mathbf{{S}}_{g}\bm{{u}}_{g}(t)$, and one obtains

[TABLE]

The absolute acceleration response of the structure can be obtained as

[TABLE]

Since accelerometers always measure absolute accelerations, therefore under the influence of seismic excitations the direct-feedthrough term $\operatorname{\mathbf{{J}}_{c}}$ vanishes.