2.1. Probabilistic model

A switching GPFLM is here introduced to enable the identification of complex nonlinearities characterized by discontinuities and/or generating different motion regimes in the dynamic response. This latent force model can be formulated as an SLDS in the following form:

$$ \left\{\begin{array}{l}p\left({\mathbf{x}}_t|{\mathbf{x}}_{t-1},{u}_{t-1},{s}_t\right)=\mathcal{N}\left({\mathbf{x}}_t|A\left({s}_t\right){\mathbf{x}}_{t-1}+B\left({s}_t\right){u}_{t-1},Q\left({s}_t\right)\right)\\ {}p\left({\mathbf{y}}_t|{x}_t,{u}_t,{s}_t\right)=\mathcal{N}\left({\mathbf{y}}_t|C\left({s}_t\right){\mathbf{x}}_t+D\left({s}_t\right){u}_t,R\left({s}_t\right)\right)\end{array}\right.{\displaystyle \begin{array}{l}\left(3\mathrm{a}\right)\\ {}\left(3\mathrm{b}\right)\end{array}} $$

where $ {s}_t $ is a switch variable denoting the active latent force model at the time instant $ t $ . For each model, the above equations correspond to the transition and observation models of the augmented linear Gaussian state-space model (LGSSM) obtained by coupling the state-space representations of the system and the latent force, as described in detail in Section 2.2.

This formulation is achieved in the same spirit of that proposed by Hartikainen and Särkkä (Reference Hartikainen and Särkkä2012) to identify dynamical systems where different latent forces can act as an input in different, and possibly overlapping, time intervals. The proposed switching latent force model mostly differs from Hartikainen and Särkkä’s model for the presence of a known input term, which enables the use of the GP latent forces for modeling the unknown nonlinear forces acting in the system. In the switching GPLFM framework, different dynamical models can be introduced to describe the evolution of the latent states and nonlinear force at different time steps, allowing the use of different governing equations to deal with different motion regimes. In the example of a dry frictional oscillator, the proposed approach allows the definition of different equations to describe the behavior of the system during sliding and sticking phases. However, to enable the presence of discontinuities in the latent nonlinear force, either at the transition between different regimes or between two sliding phases with velocities of opposite sign, it is also necessary to include a resetting model in the SLDS formulation. This model, as suggested by Hartikainen and Särkkä (Reference Hartikainen and Särkkä2012), resets the latent force components to a zero-mean Gaussian prior with a suitably chosen covariance, while leaving the output states of the system unaltered.

In the switching GPLFM formulation, the SLDS from Equation (3) is generally coupled with a discrete-time Markov model of form $ p\left({s}_t|{s}_{t-1}\right) $ , so that probabilistic inference can also be performed on the transitions between the different latent force models. If $ S $ different models are considered, with the $ S $ th model operating as a resetting model, this Markov model can be written as follows:

(4) $$ p\left({s}_t|{s}_{t-1}\right)=\left\{\begin{array}{ll}\rho, & \mathrm{if}\hskip0.35em {s}_t={s}_{t-1}\ne S\\ {}1-\rho, & \mathrm{if}\hskip0.35em {s}_t=S\hskip0.35em \mathrm{and}\hskip0.35em {s}_{t-1}\ne S\\ {}\frac{1}{S-1},& \mathrm{if}\hskip0.35em {s}_t\ne S\hskip0.35em \mathrm{and}\hskip0.35em {s}_{t-1}=S\\ {}0,& \mathrm{otherwise}\end{array}\right. $$

with $ 0<\rho <1 $ . Equation (4) can be interpreted as follows. If a certain (non-resetting) model is active at the time $ t-1 $ , there is a probability equal to $ \rho $ that such a model will remain active at the subsequent time instant $ t $ ; otherwise, a transition to the resetting model will occur. Furthermore, if the resetting model is active at the instant $ t-1 $ , in the following time step there will be a transition to any of the other models, with equal probabilities. The value of $ \rho $ must be chosen by the user, usually between 0.9 and 1 (see, e.g., Barber (Reference Barber2006); Hartikainen and Särkkä (Reference Hartikainen and Särkkä2012)).

In summary, the switching GPLFM can take into account three different sources of uncertainty in the nonlinear system identification process:

• • Nonlinear force model uncertainty. The lack of knowledge of the functional form of the nonlinear term in Equation (1) is addressed by modeling, for each motion regime, the unknown nonlinear force as a zero-mean temporal GP with a stationary kernel, which is then updated in a Bayesian manner by using the available response measurements and the relationship between nonlinear force and response determined by the physics-based model. The choice of the number of latent force models to be included in the SLDS, their mathematical formulation and the selection of appropriate kernel functions are here based on the available physical knowledge of the nonlinear system. For instance, in the case of a dry friction oscillator, two different latent force models can be used to model the dynamic behavior of the system in sliding and sticking motion regimes, respectively, and a further resetting model can be included in the SLDS to enable the presence of discontinuities in the identified nonlinear force; this is further detailed in Section 3.2. If no prior information is available regarding the type of nonlinearity, a standard GPLFM can used for a preliminary investigation of the nonlinear term. In fact, the identified latent force could give indications regarding the presence of multiple motion regimes and/or discontinuities, as in the example shown in Figure 1.

• • Measurement uncertainty. As in standard GPLFMs, the measurement noise is here assumed to follow a Gaussian distribution of covariance $ R $ ; different covariance matrices can be considered across the different latent force models. While here a zero-mean Gaussian distribution is assumed, the presence of a measurement bias could also be accounted for in the procedure; further details can be found in Barber (Reference Barber2006). In the proposed approach, the measurement noise covariance is inferred along with GP hyperparameters (see Section 2.4).

• • Regime transitions uncertainty. One of the main advantages of this approach is its ability to infer the unknown switches between the different latent force models, and therefore motion regime transitions and discontinuities. The prior knowledge about regime transitions is modeled by the Markov model introduced in Equation (4) and updated by using the response measurements, along with the latent states and forces.

The implementation of the proposed approach is schematically shown in Figure 2. In order to build the SLDS from Equation (3), it is firstly necessary to derive the state-space formulations of the system and GP latent forces, according to the procedure introduced in Section 2.2. Inference can finally be performed on the so-defined probabilistic model to retrieve the optimal GP hyperparameters and the posterior distributions for the latent states and nonlinear force, along with the sequence of motion regimes and discontinuities occurring in the time evolution. This is explained in Section 2.3.

Figure 2. Schematic representation of the switching GPLFM.

2.2. State-space formulation

Let us assume that, in Equation (1), the nonlinear term is expressed by a set of $ S-1 $ latent forces. According to the switching GPLFM formulation, one or more of these forces can be active at each time step. Each latent force can be formulated as a zero-mean temporal GP with a stationary kernel, as indicated in Equation (2). As demonstrated by Hartikainen and Särkkä (Reference Hartikainen and Särkkä2010), these GPs can be expressed as LGSSMs. Referring for simplicity to a single latent force $ f $ , this state-space representation can be expressed as follows:

(5) $$ \dot{\mathbf{f}}(t)={A}_{c,f}\mathbf{f}(t)+{L}_{c,f}w(t) $$

where the vector $ \mathbf{f}=\left[f,\overset{\;}{f},\dots, {f}^{\left(\beta -1\right)}\right] $ is constructed by taking the latent force and its derivatives, and:

(6) $$ {A}_{c,f}=\left[\begin{array}{llll}0& 1& & \\ {}& \ddots & \ddots & \\ {}& & 0& 1\\ {}-{a}_0& \dots & -{a}_{\beta -2}& -{a}_{\beta -1}\end{array}\right],\hskip0.35em {L}_{c,f}=\left[\begin{array}{l}0\\ {}\vdots \\ {}0\\ {}1\end{array}\right] $$

Following the procedure described in Hartikainen and Särkkä (Reference Hartikainen and Särkkä2010), the coefficients $ {a}_0,\dots, {a}_{\beta -1} $ , the spectral density $ q $ of the white noise process $ w(t) $ and the dimensionality $ \beta $ of the vector $ \mathbf{f} $ can be assigned so that the latent force $ f $ corresponds to the prior of a GP with the desired covariance function. The Matérn covariance functions are particularly suitable for this implementation, since they have an analytical state-space representation. The general form of this class of functions can be written as follows:

(7) $$ \kappa \left(t,{t}^{\prime}\right)={\sigma}_f^2\;\frac{2^{1-\nu }}{\Gamma \left(\nu \right)}{\left(\frac{\sqrt{2\nu }}{l}|t-{t}^{\prime }|\right)}^{\nu }{\mathcal{K}}_{\nu}\left(\frac{\sqrt{2\nu }}{l}|t-{t}^{\prime }|\right) $$

where $ {\mathcal{K}}_{\nu } $ and $ \Gamma $ are the modified Bessel and the Gamma functions respectively. $ {\sigma}_f $ and $ l $ are two hyperparameters, the first controlling the amplitude of the function, the second the lengthscale of the time variability. The smoothness parameter $ \nu $ can be selected to choose the degree of smoothness of the covariance function, ranging from the exponential ( $ \nu =1/2 $ ) to the squared exponential () $ \nu \to \mathrm{\infty} $ covariances. State-space representations have also been derived for different covariance functions, including periodic and quasi-periodic kernels (Solin and Särkkä, Reference Solin and Särkkä2014) and even some non-stationary kernels (Benavoli and Zaffalon, Reference Benavoli and Zaffalon2016). In this paper, without loss of generality, the exponential covariance function (corresponding to Matérn 1/2) will be used. In this case, Equation (5) reduces to:

(8) $$ \dot{u}(t)=-\frac{1}{l}u(t)+w(t) $$

where the white noise $ w $ is such that its spectral density is equal to $ 2{\sigma}_f^2/l $ .

The state-space representation of the GPFLM can be obtained, as suggested by Hartikainen and Särkkä (Reference Hartikainen and Särkkä2012), by coupling the state-space formulation of the zero-mean temporal GP, provided in Equation (5), with that of the dynamical system, which can be derived from Equation (1) as

(9) $$ \dot{\mathbf{z}}(t)={A}_{c,s}\mathbf{z}(t)+{B}_{c,s}\hskip0.35em \left(u(t)-f\Big(z,\dot{z}\Big)\right)=\left[\begin{array}{ccc}0& 1& \\ {}-\frac{k}{m}& -\frac{c}{m}& \end{array}\right]\left[\begin{array}{c}z\\ {}\dot{z}\end{array}\right]+\left[\begin{array}{c}0\\ {}\frac{1}{m}\end{array}\right]\left(u(t)-f\Big(z,\dot{z}\Big)\right) $$

Introducing the augmented state vector $ \mathbf{x}={\left[\begin{array}{cc}{\mathbf{z}}^{\top }& {\mathbf{f}}^{\top}\end{array}\right]}^{\top } $ , the coupling of Equations (5) and (9) results in the augmented state-space model:

(10) $$ \dot{\mathbf{x}}(t)={A}_c\mathbf{x}(t)+{B}_cu(t)+{L}_cw(t) $$

where:

(11) $$ {A}_c=\left[\begin{array}{cc}{A}_{c,s}& {B}_{c,f}\\ {}\mathbf{0}& {A}_{c,f}\end{array}\right],\hskip0.35em {B}_c=\left[\begin{array}{c}{B}_{c,s}\\ {}\mathbf{0}\end{array}\right],\hskip0.35em {L}_c=\left[\begin{array}{c}\mathbf{0}\\ {}{L}_{c,f}\end{array}\right] $$

In particular, the matrix $ {B}_{c,f} $ is obtained by adding $ \beta -1 $ columns of zeros to the matrix $ -{B}_{c,s} $ . The above stochastic differential equation (SDE) can finally be converted to a discrete-time dynamic model in the form:

(12) $$ {\mathbf{x}}_t=A{\mathbf{x}}_{t-1}+{Bu}_{t-1}+{\mathbf{w}}_{t-1},\hskip2.5em {\mathbf{w}}_{t-1}\sim \mathcal{N}(\mathbf{0},Q) $$

where the transition matrix is obtained as $ A={\exp}_m\left({A}_c\Delta t\right) $ , being $ \Delta t $ the fixed time step, and $ B={A}^{-1}\left({A}_c-I\right){B}_c $ . The covariance matrix of the process noise in the discrete model can be written as

(13) $$ Q={\int}_0^{\Delta t}{\exp}_m\left({A}_c\left(\Delta t-\tau \right)\right){L}_c{qL}_c^{\top }{\exp}_m{\left({A}_c\left(\Delta t-\tau \right)\right)}^{\top}\mathrm{d}\tau $$

where $ q $ the spectral density of the process noise in the continuous model. A more practical implementation of equation (13) is

(14) $$ Q={P}_{\infty }-{AP}_{\infty }{A}^{\top } $$

where $ {P}_{\infty } $ is the steady-state solution to the SDE describing the time evolution of the covariance of the state vector $ \mathbf{z} $ and can be retrieved by solving the Lyapunov equation (Särkkä and Solin, Reference Särkkä and Solin2019):

(15) $$ {A}_c{P}_{\infty }+{P}_{\infty }{A}_c^{\top }+{L}_c{qL}_c^{\top }=0 $$

Finally, the transition model from Equation (12) is coupled with the observation model:

(16) $$ {\mathbf{y}}_t=C{\mathbf{x}}_t+{Du}_t+{\mathbf{v}}_t,\hskip2.5em {\mathbf{v}}_t\sim \mathcal{N}(\mathbf{0},R) $$

to recover a full LGSSM, In the above equation, $ \mathbf{y} $ is the observation vector, $ C $ and $ D $ the observation matrices and $ R $ the measurement noise covariance matrix. The choice of observation matrices depends on the typology of sensor used to perform the measurements: if displacements are measured the matrices will be $ C=\left[\begin{array}{ccccc}1& 0& 0& \dots & 0\end{array}\right] $ and $ D=0 $ , while for accelerations that would be $ C=\left[\begin{array}{cccccc}-k/m& -c/m& -1/m& 0& \dots & 0\end{array}\right] $ and $ D=1/m $ . If only a single quantity is measured, the matrix $ R $ will coincide with the variance $ {\sigma}_n^2 $ of its measurement noise. Differently, if more typologies of sensors are used, the main diagonal of $ R $ will include the variances of each measurement noise, while the off-diagonal elements will represent the measurement noise covariances and will be equal to zero if noise sources are uncorrelated. The elements of the matrix $ R $ can either be inferred along with the GP hyperparameters or estimated by using different techniques (see, e.g., Bilbao et al. (Reference Bilbao, Lourens, Schulze and Ziegler2022)). The first approach will used in this paper.

The LGSSM in Equations (12) and (16) has been formulated, for simplicity, by referring to the standard GPLFM case, which is retrieved from the more general formulation from Equation (3) for $ S=1 $ . The procedure for deriving the state-space model of the GPLFM is substantially the same when more latent force models are considered. In fact, in the switching GPLFM, all the matrices of the LGSSM can be derived as described above for each latent force model $ {s}_t $ . In general, since the state-space formulations of the different latent force models are obtained independently of each other, there is no requirement for the latent forces to share the same dimensionality across the different models, implying that different kernel functions can be used in the GP priors. Moreover, not only the latent force, but also the governing equations of the system can differ for different values of $ {s}_t $ ; this will be shown in more detail in Section 3.2. Finally, the implementation of the resetting model ( $ {s}_t=S $ ), necessary to introduce discontinuities in the latent force, is achieved by formulating the transition matrix as $ A(S)=\mathrm{blkdiag}\left({\exp}_m\left({A}_{c,s}\Delta t\right),\mathbf{0}\right) $ and the covariance matrix of the process noise as $ Q(S)=\mathrm{blkdiag}\left({P}_{\infty }-{A}_s{P}_{\infty }{A}_s^{\top },{P}_0\right) $ , being $ {P}_0 $ the prior covariance selected for the latent force. This formulation, proposed by Hartikainen and Särkkä (Reference Hartikainen and Särkkä2012), enables the resetting of the latent force to its prior distribution when a discontinuity occurs without introducing discontinuities in the time evolution of displacement and velocity. It is worth underlying that, if a SDOF system is considered, only a single latent force can be included in each latent force model. The case of MDOF systems, where multiple, and possibly correlated, latent forces can be considered within the same model $ {s}_t $ , would require more complex approaches, which are beyond the scope of this contribution.

2.3. Inference of latent states and nonlinear force

Exact inference in the SLDS formulated in Equations (3) and (4) is not computationally tractable since its complexity scales exponentially with time (Bar-Shalom et al., Reference Bar-Shalom, Li and Kirubarajan2001; Lerner, Reference Lerner2002). Therefore, the filtered and smoothed posterior distributions can only be inferred in an approximated form.

In this paper, approximated inference will be performed by implementing the approach proposed by Barber (Reference Barber2006). In this approach, the filtered and the smoothed distributions of the state of each model $ {s}_t $ , are approximated by Gaussian mixtures of $ I $ and $ J $ components, respectively. The number of mixture components can be chosen by the user according to the computational budget and will generally be significantly smaller than $ {S}^t $ , which is the number of Gaussian components of the exact posterior distribution at the time instant $ t $ .

In Barber’s procedure, the filtered distributions are computed via assumed density filtering (ADF) (Alspach and Sorenson, Reference Alspach and Sorenson1972). The result of ADF at the time step $ t $ is the Gaussian mixture:

(17) $$ p\left({\mathbf{x}}_t|{s}_t,{\mathbf{y}}_{1:t},{u}_{1:t}\right)\cong \sum \limits_{i_t=1}^I{\mathcal{W}}_t\left({i}_t,{s}_t\right)\cdot \mathcal{N}\left({\mathbf{x}}_t|{\mu}_{t\mid t}\left({i}_t,{s}_t\right),{P}_{t\mid t}\left({i}_t,{s}_t\right)\right) $$

for the states of each model $ {s}_t $ and an approximated expression for the corresponding model probability $ p\left({s}_t|{\mathbf{y}}_{1:t},{u}_{1:t}\right) $ . In the above expression, the weights $ {\mathcal{W}}_t $ are obtained as $ {\mathcal{W}}_t\left({i}_t,{s}_t\right)=p\left({i}_t|{s}_t,{y}_{1:t},{u}_{1:t}\right) $ . Each step of ADF requires running $ {IS}^2 $ Kalman filters. In a similar fashion, the smoothed distributions are obtained by implementing an expectation correction (EC) algorithm, resulting in the Gaussian mixture:

(18) $$ p\left({\mathbf{x}}_t|{s}_t,{\mathbf{y}}_{1:T},{u}_{1:T}\right)\cong \sum \limits_{j=1}^J{\tilde{\mathcal{W}}}_t\left({j}_t,{s}_t\right)\cdot \mathcal{N}\left({\mathbf{x}}_t|{\mu}_{t\mid T}\left({j}_t,{s}_t\right),{P}_{t\mid T}\left({j}_t,{s}_t\right)\right) $$

where $ {\tilde{\mathcal{W}}}_t=p\left({j}_t|{s}_t,{y}_{1:t},{u}_{1:t}\right) $ , and an approximated expression for $ p\left({s}_t|{\mathbf{y}}_{1:T},{u}_{1:T}\right) $ . The EC algorithm performs $ {IJS}^2 $ RTS smoothers at each time step, and is therefore more computationally expensive then ADF. Nonetheless, it is worth mentioning that the evaluation of the marginal likelihood $ p\left({\mathbf{y}}_{1:T}\right) $ , which is required for estimating the optimal hyperparameters of the latent functions, is performed within the ADF algorithm. A detailed description of ADF and EC implementation can be found in Barber (Reference Barber2006).

2.4. Inference of GP hyperparameters

The performances of switching latent force models are not only affected by the choice of appropriate kernels for the GP latent forces, but also from the selection of the system parameters and GP hyperparameters. While an estimation procedure for uncertain system parameters will be introduced in Section 3.6, the problem of determining the optimal values for the hyperparameters in GPLFMs is addressed in what follows.

A common approach for selecting the optimal hyperparameters of the kernel function is to use maximum likelihood estimation (see, e.g., Ghahramani and Hinton (Reference Ghahramani and Hinton1996); Nayek et al. (Reference Nayek, Chakraborty and Narasimhan2019)). However, the algorithms used for determining the maximum of the likelihood function are not guaranteed to find the global maximum (Petersen et al., Reference Petersen, Oiseth and Lourens2022) and do not give information about the posterior distribution of the parameters. For this reason, different approaches have been proposed over the years. Zou et al. (Reference Zou, Lourens and Cicirello2022) and Bilbao et al. (Reference Bilbao, Lourens, Schulze and Ziegler2022) determine the optimal hyperparameters by minimizing the Hellinger distance between the empirical distribution of the measurements and the modelled Gaussian prior on the observed states; nonetheless, this approach can only be used when the response distribution is well-approximated by a Gaussian. Finally, Rogers at al. (Reference Rogers, Worden, Manson, Tygesen and Cross2018), Rogers et al. (Reference Rogers, Worden and Cross2020a), Rogers and Friis (Reference Rogers and Friis2022) recover the full posterior distribution of the parameters by using Markov Chain Monte Carlo (MCMC). In fact, while computationally more expensive, MCMC offers the advantage of a guaranteed convergence to the true posterior distribution (Gelman et al., Reference Gelman, Carlin, Rubin, Vehtari, Dunson and Stern2013).

In this paper, hyperparameters inference will be performed by using the Variational Bayes Monte Carlo (VBMC) method developed by Acerbi (Reference Acerbi2018, Reference Acerbi2020). VBMC combines variational inference and active-sampling Bayesian quadrature to perform approximate Bayesian inference, leading to a reduced computational cost compared to MCMC. Specifically, VBMC simultaneously computes an approximate posterior distribution of the model parameters and the estimated lower bound (ELBO) of the marginal likelihood. The reader is referred to the papers by Acerbi (Reference Acerbi2018, Reference Acerbi2020) for a detailed description of the VBMC method. In VBMC, the evaluation of the ELBO requires an iterative computation of the marginal likelihood of the SLDS, which is obtained as a by-product of ADF. Therefore, as schematically shown in Figure 2, each iteration requires the implementation of the LGSSM matrices by using the current values of the hyperparameters and the application of ADF. The optimal hyperparameters are finally obtained as the mean of the posterior distribution computed by VBMC and used for inferring the latent state and nonlinear force via ADF and EC.

3.1. Simulated system

Let us consider the SDOF system of mass $ m=1 $ kg, viscous damping coefficient $ c=5 $ Nsm $ {}^{-1} $ and stiffness $ k=500 $ Nm $ {}^{-1} $ shown in Figure 3a. A friction contact occurs between the mass and a ground-fixed wall, governed by a rate-dependent friction law, and a driving force $ u(t) $ is directly applied to the mass. The governing equation of this system can be formulated as:

(19) $$ m\ddot{z}+c\dot{z}+ kz+{F}_f\left(z,\dot{z}\right)=u(t) $$

Figure 3. Mass-spring-dashpot system with a friction contact between mass and ground-fixed wall (a) and its simulated response (in red) to a random phase multisine excitation (in black, divided by the stiffness) (b).

In particular, a steady-state version of the Dieterich–Ruina’s law (Dieterich, Reference Dieterich1979; Rice and Ben-Zion, Reference Rice and Ben-Zion1996), in the formulation proposed by Cabboi et al. (Reference Cabboi, Marino and Cicirello2022):

(20) $$ {F}_f\left(z,\dot{z}\right)=\left\{\begin{array}{cc}\left[{F}_{\ast }+a\ln \left(\frac{\mid \dot{z}\mid +\varepsilon }{V_{\ast }}\right)+b\ln \left(c+\frac{V_{\ast }}{\mid \dot{z}\mid +\varepsilon}\right)\right]\operatorname{sgn}\left(\dot{z}\right),& \mathrm{if}\;\dot{z}\ne 0\\ {}\hskip-19em u(t)- kz,\hskip5.879996em & \mathrm{otherwise}\end{array}\right. $$

has been used to model the friction force. The values selected for the parameters of the Dieterich–Ruina’s law are reported in Table 1. The driving force $ u(t) $ is a random phase multisine whose amplitude is generated from the JONSWAP spectrum (Hasselmann et al., Reference Hasselmann, Barnett, Bouws, Carlson, Cartwright, Enke, Ewing, Gienapp, Hasselmann, Kruseman, Meerburg, Müller, Olbers, Richter, Sell and Walden1973; Vazirizade, Reference Vazirizade2019):

(21) $$ S\left(\omega \right)=320{\left(\frac{H_s}{T_p^2}\right)}^2\cdot \frac{1}{\omega^5}\exp \left(-1.25{\left(\frac{\omega_p}{\omega}\right)}^4\right)\cdot {3.3}^{\exp \hskip0.1em \left(-\frac{1}{2{\sigma}_p^2}{\left(1-\frac{\omega }{\omega_p}\right)}^2\right)} $$

in the frequency range $ \omega =\left(0.02:0.02:100\right) $ Hz, using the parameters indicated in Table 2. The resulting input force is depicted in Figure 3b.

Table 1. Selected parameters for the steady-state Dieterich-Ruina’s law.

Table 2. Selected parameters for the JONSWAP spectrum.

It is assumed that noisy measurements of the mass displacement are available to be used as experimental observations in the GPLFM. To obtain these synthetic measurements, the dynamic response of the system has been simulated by using the event-driven variable-step Runge–Kutta (4,5) algorithm developed by Marino and Cicirello (Reference Marino and Cicirello2022), between the initial time instant $ t=0 $ and the specified final time $ t={t}_f $ . In the simulation, the initial conditions for mass position and velocity have been set to zero. The resulting displacements, shown in Figure 3b, have then been resampled by using linear interpolation according to a fixed-step time vector, whose time step is obtained from the selected sampling frequency $ {f}_s $ as $ \Delta t=1/{f}_s $ . Finally, the resulting displacement values have been polluted with artificial white noise according to a selected signal-to-noise ratio (SNR).

In the remaining of this section, the performances of the proposed switching GPLFM approach are firstly investigated and compared to the standard latent restoring force model in Section 3.2, assuming a specified set of $ {t}_f $ , $ {f}_s $ and SNR values and known physical parameters of the system. A parameter estimation approach is then proposed in Section 3.3, while the performances of the switching GPLFM for varying noise levels, observation times, and sampling frequencies are investigated in Section 3.4.

3.2. Probabilistic model for dry friction oscillators

In order to apply the proposed switching GPLFM approach to the above numerical case-study, a probabilistic model has been implemented by considering three different latent force models ( $ S=3 $ ).

• • A first model is introduced to describe the dynamic evolution of the system during the sliding phases. The implementation of this model is carried out as described in detail in Section 2.2, with the sliding friction force modeled as a zero-mean GP with an exponential (Matérn 1/2) kernel function, In fact, while the suitability of Matérn kernels for state-space implementation has been already discussed in Section 2.2, it has been verified that similar predictions are obtained if different values of the smoothness parameter (such as 3/2 and 5/2) are considered. Hence, the exponential covariance function has been selected to reduce the dimensionality of the augmented state-space model and, consequently, the computational burden.

• • The second model aims at describing the system behavior during the sticking phases. In this case, the response of the system is simply characterized by constant displacement and zero velocity, while the friction force is related to displacement and forcing function as described by Equation (20). This physical knowledge of the sticking behavior can be used to determine the prediction density as follows:

  (22) $$ p\left({\mathbf{x}}_t|{\mathbf{y}}_{1:t-1},{u}_{t-1},{i}_{t-1},{s}_{t-1},{s}_t=2\right)=\mathcal{N}\left({\mathbf{x}}_t|{\mu}_{t\mid t-1}\left({i}_{t-1},{s}_{t-1}\right),{P}_{t\mid t-1}\left({i}_{t-1},{s}_{t-1}\right)\right) $$
  where:
  (23) $$ {\mu}_{t\mid t-1}\left({i}_{t-1},{s}_{t-1}\right)={\left[\begin{array}{c}\unicode{x1D53C}\left({z}_{t-1\mid t-1}\left({i}_{t-1},{s}_{t-1}\right)\right),0,{u}_{t-1}-k\cdot \unicode{x1D53C}\left({z}_{t-1\mid t-1}\left({i}_{t-1},{s}_{t-1}\right)\right)\end{array}\right]}^{\top } $$
  and $ {P}_{t\mid t-1}\left({i}_{t-1},{s}_{t-1}\right)={P}_{t-1\mid t-1}\left({i}_{t-1},{s}_{t-1}\right) $ . If a different kernel from Matérn 1/2 is considered, zeros can be added at the end of $ {\mu}_{t\mid t-1}\left({i}_{t-1},{s}_{t-1}\right) $ to obtain the required dimensionality. It is worth underlying that the prediction density from Equation (22) can be directly filtered and smoothed in the present form, without need of defining the matrices $ A\left({s}_t\right) $ and $ B\left({s}_t\right) $ for the sticking model.

• • The third model is the resetting model, which is implemented as explained in Section 2.2, setting the prior covariance of the latent force to $ {P}_0=0.05 $ . When activated, this model resets the latent friction force to zero, with $ {P}_0 $ as covariance), while leaving displacement and velocity unaltered. Since, as described by Equation (4), this model can only remain active for a single time step, the result is that the values assumed by the friction force immediately before and after the reset will not directly affect each other. This allows for the presence of a discontinuity in the friction force, either between a sticking and a sliding phases or between two sliding phases with different velocity sign.

Finally, the model transitions are governed by the Markov model introduced in Equation (4), where $ \rho $ has been set to 0.92.

3.3. Inference of latent states and nonlinear force with known system parameters

Nonlinear force identification has firstly been performed by considering the simulated dynamic response of the system introduced above for $ {t}_f=5\hskip0.1em $  s and $ {f}_s=500 $  Hz. The mass displacement has been corrupted with artificial white noise by using the AWGN function of Matlab (2023), setting SNR $ =80 $ . The so-produced white noise has a standard deviation of $ 7.8623\times {10}^{-3}\hskip0.1em $  mm, corresponding to a variance $ {\sigma}_n^2=6.1815\times {10}^{-5}{\mathrm{mm}}^2 $ .

In this example, the system parameters have been assumed as known a priori and set to their true values indicated at the beginning of this section. The posterior distributions of the hyperparameters $ {\sigma}_f^2 $ , $ l $ and $ {\sigma}_n^2 $ , whose assigned priors are reported in Table 3, have been inferred by using VBMC. The prior and posterior distributions of the hyperparameters are shown in Figure 4, while the convergence of the ELBO is reported in Figure 5. As specified in Section 2.3, in the iterative evaluation of the ELBO, the log marginal likelihood is obtained as by-product of the ADF. Once the optimal GP hyperparameters have been estimated as the mean values of their posterior distributions, the latent states and nonlinear restoring force have been inferred by implementing ADF and EC. The number of Gaussian mixture components in the filtering and smoothing passes has here been set to $ I=J=3 $ . The estimated displacement, velocity, friction force and regime transitions are reported, along with the ground truth, in Figures 6 and 7. The agreement between the mean of the inferred distribution and the ground truth can be quantified by evaluating the normalized mean squared error (NMSE), which can be defined as

(24) $$ \mathrm{NMSE}\left[g\right]=\frac{1}{T{\sigma}_g^2}\sum \limits_{t=1}^T{\left({g}_t-\unicode{x1D53C}\left({\hat{g}}_t\right)\right)}^2 $$

where $ {g}_t $ and $ \unicode{x1D53C}\left({\hat{g}}_t\right) $ are the ground truth and the estimated mean value of a certain quantity $ g $ at the time step $ t $ , respectively, $ {\sigma}_g^2 $ is the variance of $ g $ and $ T $ is the number of data points. According to the above definition, a value of zero corresponds to a perfect fit between true and estimated functions, while a unitary value is obtained if $ \unicode{x1D53C}\left({\hat{g}}_t\right) $ is equal to the mean value of the ground truth at every point. The very low NMSE scores reported in Table 4 for displacement and velocity confirm their excellent agreement with the ground truth. The nonlinear force identification score is NMSE $ \left[{F}_f\right]=1.95\% $ , which is generally considered a very good result. Also, the acceleration, which is determined as a linear combination of states and friction force, presents a very good agreement with the ground truth (NMSE $ \left[\ddot{z}\right]=0.46\% $ ), despite the presence of discontinuities in its time evolution. Finally, in the bottom frame of Figure 6, it can be noted that most discontinuities and regime transitions are well captured by the model transitions in the switching GPLFM. Among the seven stops occurring in the observed time window, only the second has been missed; however, it can be noted that this sticking phase is particularly short and does not affect significantly the friction force value. Similarly, spikes corresponding to the resetting model activation can be observed every time that discontinuities occur between the sliding stages. It is worth remembering that the switching GPLFM frameworks allows for more latent force models to be active simultaneously and, therefore, the blue line in Figure 6 is only referred to the predominant model at each time step. A full picture of the model probabilities is instead provided in Figure 7. Here, it can be seen that the resetting model is also activated in the transitions between sliding and sticking phases, although its probability does not always overcome that of the sticking model.

Table 3. Prior distributions for the GP hyperparameters used in the numerical case-study.

Figure 4. Prior and posterior distributions of the GP hyperparameters inferred by VBMC in switching GPLFM $ (I=J=3) $ .

Figure 5. Convergence of the evidence lower bound (ELBO) in VBMC inference for switching GPLFM $ (I=J=3) $ .

Figure 6. Latent states, acceleration, nonlinear friction force, and model sequence inferred by switching GPLFM ( $ I=J=3 $ ) vs ground truth. Models 1, 2, and 3 stand for sliding, sticking, and resetting models, respectively.

Figure 7. Model probabilities estimated by switching GPLFM ( $ I=J=3 $ ). Models 1, 2, and 3 stand for sliding, sticking, and resetting models, respectively.

Table 4. Latent states, acceleration and nonlinear force identification error scores, and optimal hyperparameters for standard and switching GPLFMs with a varying number of Gaussian mixture components. The ground truth for the measurement noise variance is $ {\sigma}_n^2=6.182\times {10}^{-11}{m}^2 $ .

3.4. Comparison between standard and switching GPLFMs

In this subsection, the performance of the switching GPLFM in latent states and nonlinear force identification are investigated for varying number of Gaussian mixture components. For simplicity, the same number of Gaussians is used for ADF and EC $ \left(I=J\right) $ , although the EC implementation, as proposed by Barber (Reference Barber2006), also allows for a smaller number of Gaussians in the smoothed distribution ( $ J<I $ ). The standard GPLFM, which can simply be seen as a particular case of switching GPLFM for $ S=1 $ , $ I=1 $ and $ J=1 $ , will be here considered as a reference case. In this case, ADF and EC will automatically reduce to Kalman filtering and RTS smoothing, respectively.

The latent states and nonlinear force identification errors obtained by using standard and switching GPLFMs (with 3, 5, 7, and 9 Gaussian components) are reported in Table 4, along with the optimal GP hyperparameters. Among the error scores, the normalized mean variance (NMV) index has been introduced as:

(25) $$ \mathrm{NMV}\left[g\right]=\frac{1}{T{\sigma}_g^2}\sum \limits_{t=1}^T\unicode{x1D54D}\left({\hat{g}}_t\right) $$

According to this definition, a zero value is obtained if the variance of the predicted function is zero at every point, while a unitary NMV corresponds to a predictive function whose variance at every time step $ t $ is equal to the overall variance of the ground truth. The identified nonlinear forces are also shown in Figure 8, where they are graphically compared to the true friction force in the time domain; a small section of the observed time window is enlarged on the right to illustrate how the different GPLFM approaches perform at discontinuities.

Figure 8. Nonlinear friction force inferred by standard and switching GPLFMs for varying number of Gaussian components in ADF and EC ( $ I=J $ ) vs ground truth.

Figure 8a shows a good overall agreement between true and identified friction force obtained by using standard GPLFM. While discontinuities cannot be captured by a single latent force, the mean prediction follows quite closely the ground truth, resulting in an identification error of NMSE $ \left[{F}_f\right]=3.27\% $ . On the other hand, it can be observed that the confidence bounds around the mean prediction remain quite large (NMV $ \left[{F}_f\right]=2.85\% $ ). The optimal hyperparameters reported in Table 4 suggest that the good agreement between the inferred and true friction force has been obtained by reducing the length-scale hyperparameter of the kernel function. In fact, small values of $ l $ allow for faster variations of the latent function, increasing its capability to adapt to the sharp variations of the friction force. Therefore, this partially compensate the inability of a single GP latent force to model discontinuities and different motion regimes (i.e., sliding and sticking phases). Nonetheless, short length-scales are also associated with well-known downsides, such as the increase of model complexity and a more pronounced tendency to overfitting, particularly if noisy inputs are considered (Rasmussen and Williams, Reference Rasmussen and Williams2006; Beckers, Reference Beckers2021).

In the remaining frames of Figure 8, it can be observed that the introduction of the sticking and resetting models enable the capture of most discontinuities and transitions between different motion regimes in the friction force. This leads to a significant improvement of the nonlinear force identification performance: as reported in Table 4, the NMSE score associated with the friction force drops to $ 1.95\% $ when three Gaussian components are considered and even to $ 1.66\% $ when $ I=J=5 $ . Nonetheless, further increases in the number of components do not appear to lead to further improvements and NMSE has been observed to increase above $ I=J=7 $ . A possible explanation of this trend can be found in the presence of some increasingly irregular transition between motion regimes as the number of Gaussians is increased (see, for example, the transition at $ t\cong 2.4\hskip0.35em $ s). This result is also reflected by the NMV score; nonetheless, the standard deviations of the friction force are overall significantly smaller compared to the case of the standard GPLFM. Finally, looking at the optimal hyperparameters, it can be noted that the length scales obtained by introducing the switching GPLFM are significantly larger compared to the value obtained for the standard latent restoring force model, while overall good performances are observed in the measurement noise estimation independently of the number of Gaussian components.

3.5. Nonlinear friction force characterization

In the previous subsections, the application of the switching GPLFM to the simulated response of the dry friction oscillator has enabled the identification of the time evolution of the nonlinear friction force, state, and regime transitions. One of the advantages of GPLFMs is that no specific assumptions are required regarding the functional form of the latent nonlinear force. In this case-study, basic physical constraints have been considered in the implementation of the sticking regime, by imposing that stops are characterized by a constant displacement and zero velocity. However, no assumptions have been made about the dependence of the friction force on displacement and velocity in sliding conditions. Here it is shown how the identified friction force and state, along with the inferred stick–slip regime transitions, can be used to reconstruct the underlying friction law. This law, along with a correct estimation of the system parameters (dealt with in Section 3.6), enables the implementation of a robust forward model which can be used to predict the response of the dry friction oscillator to different inputs. The reconstruction of the friction model includes the determination of (i) the friction force–velocity relationship in sliding conditions and (ii) the value of the static friction force. These two steps are presented in what follows (Figure 9).

Figure 9. Nonlinear friction force vs mass displacement (a) and velocity (b): comparison between simulated (in red) and inferred (in blue) values. The fitted friction force–velocity curve is also reported (in green) in (b).

The characterization of the friction force–velocity relationship, or more in general of the sliding friction model, can simply be performed by fitting the latent force-latent states estimates, as shown in Figure 9. Several different approaches may be used for fitting; for instance, the parameters of an existing friction model can be estimated by minimizing the least squared error, or a black-box approach, such as a neural network or a further GP, can be considered. For instance, in Rogers and Friis (Reference Rogers and Friis2022), where a standard GPLFM is used for the identification of a Duffing oscillator, the nonlinear term is characterized by using a Bayesian Information Criterion to establish the most likely polynomial order of the nonlinear force–displacement relationship. Here, it has been chosen to fit the nonlinear force–velocity estimates inferred by the switching GPLFM ( $ I=J=5 $ ) with the steady-state Dieterich-Ruina’s law introduced in Equation (20a). This choice is motivated by the flexibility of this friction model in describing, for varying parameters, both the velocity-weakening and velocity-strengthening behaviors typically observed in experimental tests (Cabboi et al., Reference Cabboi, Marino and Cicirello2022); obviously, different laws or functions could be considered. The optimal values of the parameters $ a $ , $ b $ , $ c $ and $ {F}_{\ast } $ have been determined by minimizing the least squared error, while $ {V}_{\ast } $ and $ \varepsilon $ have been set to the values reported in Table 1. The resulting fitted friction force-velocity curve is illustrated in Figure 9b.

While the above fitting procedure could also be applied to estimates from the standard GPLFM, the information provided by the switching GPLFM about regime transitions also enables the estimation of the static friction force. In fact, it is well-known that the transitions from sticking to sliding regime take place when the instantaneous absolute value of the friction force value equates to the static friction force (see Cabboi et al. (Reference Cabboi, Marino and Cicirello2022) for reference). Therefore, following the model sequence illustrated in Figure 6, it is possible to use the absolute values of the latent force in correspondence of these transitions as estimates of the static friction force. To finalize the procedure, the mean and the standard deviation of these estimates are evaluated, as graphically shown in Figure 10. The estimated static friction force $ {\hat{F}}_s $ can be included in the fitted friction model by imposing $ {F}_f\left(z,0\right)={\hat{F}}_s $ in Equation (20a). This position results in the constraint:

(26) $$ b=a+\frac{{\hat{F}}_s-{F}_{\ast }}{\ln \left({V}_{\ast}\right)-\ln \left(\varepsilon \right)} $$

leaving only $ a $ , $ c $ and $ {F}_{\ast } $ as independent parameters of the fitting model. Figure 10 illustrates the comparison of the estimated static friction force, reported along with its confidence bounds, and the true value selected in the numerical simulation. It can be observed how those values, reported for $ I=J=5 $ are very close, with the ground truth falling well within the bounds. It is worth underlining that the accuracy of this procedure strongly depends on the observed number of stops; in a continuously sliding motion, it will not be possible to make any predictions on the static friction value.

Figure 10. Estimation of the static friction force from the switching GPLFM results.

3.6. Parameter estimation

In the previous analyses, the mass, the viscous damping coefficient and the spring stiffness of the system have been assumed as known. However, it is often the case in experimental contexts that the values of these physical parameters are not exactly known. In several applications of the GPLFM to mechanical systems (see, e.g., Rogers et al., Reference Rogers, Worden and Cross2020a, Reference Rogers, Worden and Cross2020b), the uncertain system parameters are inferred, along with the GP hyperparameters, by using sampling approaches such as MCMC. This approach also requires to provide a prior distribution for these parameters, which is usually specified based on engineering and physical knowledge. Nonetheless, as highlighted by Rogers and Friis (Reference Rogers and Friis2022), the presence of a latent nonlinear force can introduce a bias in the parameter estimation. For instance, the presence of a Duffing-type nonlinearity is likely to lead to incorrect linear stiffness estimates. In return, incorrect parameter estimates will also alter the identified nonlinear force.

In the case of a dry friction oscillator, it has been observed that parameter inference, performed by using VBMC, does not lead to accurate estimates when either standard or switching GPLFM is applied, even if informative prior distributions are provided. Therefore, it can be concluded that, in this case, there is no real advantage in performing parameter inference rather than simply assuming initial guesses based on engineering insights. A procedure for correcting these initial guesses, based on the inferred latent force and the physical knowledge of the system, is presented in what follows.

Let us denote with $ m $ , $ c $ , $ k $ the true values of the system parameters, and with $ \hat{m} $ , $ \hat{c} $ , $ \hat{k} $ the initial guesses, so that:

(27) $$ m=\hat{m}+\Delta m\hskip3em c=\hat{c}+\Delta c\hskip3em k=\hat{k}+\Delta k $$

By substituting Equation (27) into Equation (19), it can be seen how using incorrect parameter values will introduce additional linear forces in the governing equation of the system:

(28) $$ \hat{m}\ddot{z}+\hat{c}\dot{z}+\hat{k}z+\underset{{\mathrm{F}}_{\mathrm{L}}}{\underbrace{F_f+\Delta m\cdot \ddot{z}+\Delta c\cdot \dot{z}+\Delta k\cdot z}}=u(t) $$

Therefore, if the latent states $ z $ , $ \dot{z} $ and the latent force $ {F}_L $ are inferred by considering incorrect initial guesses as system parameters, these additional forces will become part of the latent force, alongside the nonlinear friction force. The acceleration $ \ddot{z} $ , which is not directly inferred when applying the GPLFM, can be retrieved, from the above equation, as

(29) $$ \ddot{z}=\frac{1}{\hat{m}}\left(u(t)-\hat{k}z-\hat{c}\dot{z}-{F}_L\right) $$

Substituting Equation (29) into the expression of the latent force $ {F}_L $ and rearranging, it is obtained that:

(30) $$ {F}_L=\left(\frac{\hat{m}}{\hat{m}+\varDelta m}\right){F}_f+\left(\frac{\hat{m}}{\hat{m}+\varDelta m}\varDelta k-\frac{\varDelta m}{\hat{m}+\varDelta m}\hat{k}\right)z+\left(\frac{\hat{m}}{\hat{m}+\varDelta m}\varDelta c-\frac{\varDelta m}{\hat{m}+\varDelta m}\hat{c}\right)\dot{z}+\left(\frac{\varDelta m}{\hat{m}+\varDelta m}\right)u $$

The above equation highlights how incorrect physical parameters introduce a linear dependency of the inferred latent force on the inferred states and known forcing function. In particular, it is worth noting that, if a correct estimate is provided for the mass ( $ \Delta m=0 $ ), the viscous damping and stiffness errors would simply introduce an additional linear trend in the latent force–velocity and latent force–displacement, respectively. Differently, the presence of a mass error term not only alters these linear trends, but also introduce a dependency on the forcing function and a scaling effect of the friction force.

The latent states and force have been inferred for the dry friction oscillator by considering the incorrect parameters $ \hat{m} $ , $ \hat{c} $ and $ \hat{k} $ reported in Table 5. The relationships between the latent force estimates and the displacement, velocity and forcing function are illustrated in Figure 11, respectively. In these figures, the presence of the aforementioned linear trends can be observed, along with a general alteration of the latent force patterns due to the presence of a significant mass error. To quantify the linear dependencies expressed in Equation (30), the latent force can be fitted with a multidimensional linear function, whose coefficients can be retrieved by using a least squared error approach. In order to use all the available data, the symmetry of the friction force $ {F}_f $ , and consequently of the latent force $ {F}_L $ , with respect to the origin ( $ z=0 $ , $ \dot{z}=0 $ , $ u=0 $ ) can be exploited for superposing its positive and negative branches, so that they can finally be fitted with the linear function:

(31) $$ {F}_L={A}_0+{A}_1z+{A}_2\dot{z}+{A}_3u $$

Table 5. True, guessed and corrected values of the physical parameters $ m $ , $ c $ and $ k $ of the dry friction oscillator. The relative errors of guessed and corrected parameters are referred to the true value.

Figure 11. Nonlinear friction force vs displacement (a), velocity (b), and driving force (c): comparison between ground truth (in red) and values inferred by using incorrect physical parameters (in blue).

By comparing Equations (30) and (31), it is possible to relate the parameter errors to the evaluated linear coefficients as

(32) $$ \Delta k=\frac{A_1+{A}_3\hat{k}}{1-{A}_3}\hskip3em \Delta c=\frac{A_2+{A}_3\hat{c}}{1-{A}_3}\hskip3em \Delta m=\frac{A_3\hat{m}}{1-{A}_3} $$

The estimated parameter errors can finally be added to the initial guesses as indicated in Equation (27) to retrieve the true parameters of the system. The corrected parameter estimates for the dry friction oscillator case-study are reported in Table 5. At this stage, the friction force can be retrieved from the latent force as

(33) $$ {F}_f=\frac{1}{1-{A}_3}\left({F}_L-{A}_1z-{A}_2\dot{z}-{A}_3u\right) $$

However, it is worth noting that, after the correction of the system parameters, the optimal GP hyperparameters might be different from those initially estimated from VBMC. Therefore, depending on the available computational budget and the relevance of the parameter corrections, it might be preferable to newly infer the optimal hyperparameters and thus the latent states and nonlinear force. This latter approach has been followed in this paper. The inferred friction force is plotted in Figure 12 as a function of the displacement, velocity, and driving force, respectively, and presents a very good agreement with the ground truth. As reported in Table 5, the residual relative errors between the corrected and true parameters are all negligible, despite the very large initial errors. The viscous damping coefficient is the only parameter whose residual error is above 1%. Although this error would be acceptable for most applications, it is worth mentioning that this slight overestimation is related to the assumed rate-dependent friction model. In fact, as shown in Figure 12b, the selected Dieterich-Ruina’s law exhibits an increasing trend with the sliding velocity, which has been identified as an additional viscous damping term.

Figure 12. Nonlinear friction force vs displacement (a), velocity (b), and driving force (c): comparison between ground truth (in red) and values inferred by using the corrected physical parameters (in blue).

3.7. Performance analysis for varying noise levels, observation times, and sampling frequencies

In the previous subsections, the nonlinear system identification performances of the switching GPLFM application to a dry friction oscillator have been investigated by considering a set of observations with a specific level of measurement noise and a fixed number of samples. Therefore, a further step of this investigation consists in analyzing how the proposed approach performs when different noise levels, signal durations, or sampling frequencies are taken into account. The switching GPLFM is thus applied to the case study investigated in this section for varying $ {\sigma}_n $ , $ {t}_f $ and $ {f}_s $ , assuming the initial guess reported for the system parameters in Table 5 and selecting $ I=J=3 $ . The following performance estimators are considered:

• • the nonlinear force identification error is evaluated as the NMSE of the identified friction force with respect to the ground truth (see Equation (24));

• • the prediction error aims at evaluating the predictive capabilities of the forward model implemented by considering the best-fitting Dieterich’s–Ruina friction law and the corrected parameters estimates. It is calculated as the NMSE of the mass displacement obtained numerically from the forward model with respect to the ground truth;

• • the motion regime identification error is calculated as the relative error between the true and the identified motion regimes. The most probable latent force model is considered as active model at every time step;

• • the parameter estimation errors are the relative errors between estimated and true physical parameters.

The switching GPLFM performances have been investigated for seven different levels of measurement noise (SNR = [60, 70, 80, 90, 100, 125, and 150] dB) and for a noise-free measurement (SNR = $ \infty $ ), while maintaining the final time and the sampling frequency of the observations set to $ {t}_f=5\hskip0.1em $  s and $ {f}_s=500 $  Hz. The corresponding values of the standard deviation $ {\sigma}_n $ are reported in Table 6. The performance indexes are plotted in Figure 13. In this figure, it is possible to observe that switching GPLFM presents similarly good performances with respect to all the measured errors when the signal-to-noise ratio is equal or above to 80 dB. For higher noise levels, the main effect appears to be a decrease of the motion regime identification accuracy, which also leads to lower identification scores. Nonetheless, all the scores maintain generally acceptable value for SNRs larger than 60 dB.

Table 6. Conversion between signal-to-noise ratio and standard deviation values for the measurement noise.

Figure 13. Nonlinear system identification performances of the switching GPLFM ( $ I=J=3 $ ) applied to the dry friction oscillator case-study for $ {t}_f=5\hskip0.1em $ s, $ {f}_s=500 $ Hz and varying measurement noise levels.

A further performance analysis has been carried out maintaining the noise levels and the sampling frequency set to SNR  $ =80 $  dB and $ {f}_s=500 $  Hz while varying the simulation time between 2 and 20 s. The resulting error indexes, reported in Figure 14, are substantially constant across different observed times. Therefore, it can be deduced that the switching GPLFM can offer stable performances for a varying number of samples, as long as the time step among the observations is maintained constant.

Figure 14. Nonlinear system identification performances of the switching GPLFM ( $ I=J=3 $ ) applied to the dry friction oscillator case-study for SNR $ =80 $  dB, $ {f}_s=500 $  Hz, and varying simulation times.

The last investigation has been performed by varying the sampling frequency of the measurements between 100 and 2000 Hz; this corresponds to a variation of the fixed time step between 0.0005 and 0.01 s. The SNR and time span of the measurements have been set to 90 dB and 5 s, respectively. Observing the resulting scores, illustrated in Figure 15, it is clear that the sampling frequency has a larger impact on the performances than the simulation time. While large errors are displayed for most of the investigated parameters when $ {f}_s=100 $  Hz, indicating that the switching GPLFM cannot perform well at very low sampling frequencies, very good scores are obtained when the sampling frequency is increased to 250 Hz. At the other of the investigated range, it is very interesting to observe that the nonlinear force identification error increases significantly at high sampling frequency, indicating the oversampling issues might occur when applying the GPLFM. On the other hand, the prediction, motion regime identification and parameter estimation errors do not appear to be affected by oversampling.

Figure 15. Nonlinear system identification performances of the switching GPLFM ( $ I=J=3 $ ) applied to the dry friction oscillator case-study for SNR $ =80 $  dB, $ {t}_f=5\hskip0.1em $  s and varying sampling frequencies.

4.1. Test rig and mechanical model

The experimental tests have been carried out on the rig shown in Figure 16a, which is briefly described in what follows. The main structure is a single-storey frame made up of two metal plates connected by four doubly-bolted thin metal bars. The bottom plate is excited by the alternate motion of a shaft, which is connected to an electric motor through a Scotch-yoke mechanism. This results in an approximatively harmonic base excitation, whose driving frequency can be specified by selecting the rotation speed of the motor. A brass disc is placed on the top steel plate to create a friction contact. The disc is mounted on a counterweight system pinned to the external frame, so that the normal force exerted on the top plate, and therefore the friction force amplitude, can be adjusted by shifting the weights along the counterweight axis. The displacement of the two plates is measured by laser position sensors, thus providing input and output measurements during the tests. A more detailed explanation of the setup and testing procedure can be found in Marino and Cicirello (Reference Marino and Cicirello2020).

Figure 16. Picture (a) and mechanical model (b) of the test rig: base-excited single-store frame with a brass-to-steel contact.

In the above publication, it was shown that the dynamic behavior of the single-storey frame is well reproduced by the single-degree-of-freedom mass-spring model shown in Figure 16b when driving frequencies up to at least 8 Hz are considered. Moreover, the experimental results obtained by Marino and Cicirello for the dynamic response of this structure agree well with the theoretical results obtained analytically and numerically by considering the Coulomb friction model. Since Coulomb’s law introduces a discontinuity in the friction force at zero sliding velocity, it is reasonable to assume that the single-storey frame behaves as a discontinuous nonlinear system, thus representing a suitable case-study for the proposed method.

The governing equation of the mechanical system schematized in Figure 16b can be written as

(34) $$ m\ddot{z}+c\dot{z}+ kz+{F}_f\left(z,\dot{z}\right)= ku(t)+c\dot{u}(t) $$

where $ u(t) $ and $ \dot{u}(t) $ are the displacement and velocity of the base. Differently from the mechanical system investigated in Section 3, where the known driving force was directly applied to the mass, the input term of Equation (34) is made up of two terms, both dependent on generally unknown physical parameters, and the derivative of input base excitation must also be provided. The probabilistic model implemented in Section 3.2 can be easily adapted to deal with a base-excited oscillator. Specifically:

• • in the state-space model of the system, the input vector will be written as $ \mathbf{u}={\left[\begin{array}{cc}u& \dot{u}\end{array}\right]}^{\top } $ and the matrix $ {B}_{c,s} $ will be formulated as

  (35) $$ {B}_{c,s}=\left[\begin{array}{cc}0& 0\\ {}-\frac{k}{m}& -\frac{c}{m}\end{array}\right] $$

• • in the predicted mean value of the latent restoring force from Equation (23), the term $ {u}_{t-1} $ will be replaced by $ {ku}_{t-1}+c{\dot{u}}_{t-1} $ .

Nonetheless, while the formulation of the switching latent force model is not heavily affected by the presence of the base excitation, difficulties may arise regarding the estimation of the physical parameters and input vector.

The first problem is related to the non-identifiability of the mass. In fact, rewriting Equation (34) as:

(36) $$ \ddot{z}=\frac{k}{m}u(t)+\frac{c}{m}\dot{u}(t)-\frac{k}{m}z-\frac{c}{m}\dot{z}-\frac{1}{m}{F}_f\left(z,\dot{z}\right) $$

it can be observed that, in the right-hand side, all the terms except the unknown friction force are multiplied by either $ k/m $ or $ c/m $ . Therefore, if the mass value is changed from its true value $ m $ to an incorrect estimate $ \hat{m} $ while keeping these ratios constant, the only effect on the identified states and latent force is that the latter would be scaled by a factor $ \hat{m}/m $ . Obviously, since the friction force amplitude is unknown, this scaling factor cannot be determined from the identified latent force. In conclusion, independently of the system identification approach considered, it is not possible to identify the mass parameter of a base-excited system from its base and mass displacements. In this contribution, the value $ m=3.0799 $  kg will be assigned to the mass. This value has been obtained by performing a hammer test on the top plate of the frame after removing the counterweight system and, therefore, the friction contact. The experimental frequency response function has thus been fitted with that of a mass-spring-damper system in the frequency range 0–6 Hz, retrieving the optimal parameters of the system with a least squares method. The so-obtained physical and modal parameters of the system are reported in Table 7, where they are compared to the parameters estimated by fitting the identified latent force with the linear function:

(37) $$ {F}_L={A}_0+{A}_1\left(z-u\right)+{A}_2\left(\dot{z}-\dot{u}\right) $$

and correcting the initial estimates by adding $ \Delta k={A}_1 $ and $ \Delta c={A}_2 $ , according to the procedure from Section 3.3.

Table 7. Physical and modal parameters of the single-storey frame: estimated values from hammer testing vs initial guess and corrected values in the identification procedure.

A further issue in the application of the switching GPLFM to this experimental case-study is that the base displacement and velocity must be provided as an input. However, although the fundamental frequency of the base excitation is selected by the user, only noisy measurements of the bottom plate displacements are available from the laser sensors. To recover the driving functions $ u(t) $ and $ \dot{u}(t) $ , the following state-space model has been implemented from Newton’s laws of motion:

$$ \hskip12.1em \left\{\begin{array}{ll}\left[\begin{array}{c}{u}_t\\ {}{\dot{u}}_t\end{array}\right]& \hskip-.75em =\left[\begin{array}{cc}1& \Delta t\\ {}0& 1\end{array}\right]\left[\begin{array}{c}{u}_{t-1}\\ {}{\dot{u}}_{t-1}\end{array}\right]+\left[\begin{array}{c}\frac{1}{2}\Delta {t}^2\\ {}\Delta t\end{array}\right]{a}_{t-1},\hskip4em {a}_{t-1}\sim \mathcal{N}\left(0,{\sigma}_a^2\right)\hskip9.7em \left(38\mathrm{a}\right)\\ {}\hskip1.5em {y}_{u,t}& \hskip-.75em =\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}{u}_t\\ {}{\dot{u}}_t\end{array}\right]+{v}_{u,t},\hskip13.5em {v}_t\sim \mathcal{N}\left(0,{\sigma}_u^2\right)\hskip11.15em \left(38\mathrm{b}\right)\end{array}\right. $$

by assuming that an unknown constant acceleration $ {a}_t $ , modelled as a zero-mean normal distribution with variance $ {\sigma}_a^2 $ , acts between the time steps $ t $ and $ t+1 $ . The above state-space model has then been computed via Kalman filtering and RTS smoothing, by considering the optimal values of the hyperparameters inferred by VBMC. It is worth mentioning that, while the proposed procedure can be conveniently applied in the assumption that the base motion is not affected by either the dynamic response of the structure or the friction force, a more general approach would consist in coupling the above state-space model with the augmented state-space representation of the latent force models.

4.2. Inference and results

Nonlinear system identification has been performed by applying the proposed switching GPLFM ( $ I=J=3 $ ) to the experimental data obtained by setting a fundamental driving frequency of 1 Hz and a normal force equal to 5.5 N. The displacements of the two plates, reported in Figure 17, have been recorded for 10 s with a sampling frequency of 250 Hz.

Figure 17. Measured displacements of the bottom and top plates of the single-storey frame under harmonic excitation. The driving frequency is set at 1 Hz and a normal load of 5.5 N is applied to the top plate.

The estimated physical and modal parameters are listed in Table 7. In particular, the estimated stiffness and natural frequency agree well with the measured values from hammer testing, while the estimated viscous damping and the corresponding damping ratio are slightly overestimated. Nonetheless, viscous damping has negligible values in both cases, in agreement with the results presented by Marino and Cicirello (Reference Marino and Cicirello2020).

The identified latent states and friction force are reported in Figure 18, along with the sequence of model probabilities. In this figure, it is possible to observe that the estimated displacements present an excellent agreement with the experimental observations, scoring an NMSE equal to 4.889 $ \times {10}^{-6} $ . The displacement of the top plate is a quite regular two-stops stick–slip motion. As shown in the bottom frame of Figure 18, the transitions between the sliding and sticking phases are generally well captured by the switching latent force model. Nonetheless, during the stops occurring at $ z<0 $ in the top plate motion, slow variations can be observed before the onset of the subsequent sliding phases. This particular behavior clearly renders more difficult the identification of the motion regime, leading to the presence of grey areas in the model probabilities corresponding to these stops. The identified friction force does not appear to be affected by the larger uncertainty that characterizes negative stops; in fact, regular patterns are observed during each sticking phase. The sliding friction force is also characterized by generally regular behaviors, with a slight decrease in amplitude during each sliding phase. This trend can also be visualized in the latent friction force–velocity point estimates reported (with magenta dots) in Figure 19.

Figure 18. Latent states, acceleration, nonlinear friction force, and model probabilities inferred by switching GPLFM ( $ I=J=3 $ ). Models 1, 2, and 3 stand for sliding, sticking, and resetting models, respectively.

Figure 19. Friction force vs velocity for varying applied normal loads: points estimates (dots), fitted friction laws (continuous lines) and $ \pm 3\sigma $ confidence intervals (shaded areas).

The friction force–velocity estimates have been fitted with the steady-state Dieterich-Ruina’s law introduced in Equation (20a); the value of the static friction force has been determined as described in Section 3.5. In Figure 19, it can be observed that both branches of the friction force–velocity curve follow a monotonic behavior, characterized by a decrease of the friction force amplitude for increasing absolute values of the sliding velocity. This particular pattern, known in the friction literature as Stribeck effect, is reproduced by the Dieterich–Ruina’s law due to the very small identified values for the parameters $ a $ and $ c $ (see also Cabboi et al. (Reference Cabboi, Marino and Cicirello2022)).

To evaluate the predictive capabilities of the switching GPLFM, the dynamic response of the system has been simulated numerically by considering the fitted friction model and the estimated physical parameters (see Table 7). The measured and simulated mass displacements present a very good agreement, also visible in the top frame of Figure 18; their comparison yields an NMSE score of 0.4316%.

Finally, the proposed approach has been applied to experimental data obtained for different normal force levels, ranging from 2.5 to 5.5 N. The obtained friction force–velocity estimates, along with the corresponding fitted friction laws, are reported in Figure 19. It can be observed that the friction force–velocity curves present similar trends across the different applied normal forces. In addition, the curves are almost equally spaced, suggesting that the friction force amplitude increases linearly with the normal force, as expected. The predictive performances of the fitted friction models do not present significant variations are substantially maintained for the different values of $ N $ , with NMSE scores never exceeding 1% among the observed cases.