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

HAL
open science

# Identification and monitoring of stochastic linear subsystems with unknown local nonlinearities via output injection

Neha Aswal, Adrien Mélot, Laurent Mevel, Qinghua Zhang

► To cite this version:

Neha Aswal, Adrien Mélot, Laurent Mevel, Qinghua Zhang. Identification and monitoring of stochastic linear subsystems with unknown local nonlinearities via output injection. Mechanical Systems and Signal Processing, 2025, 238, pp.113239. 10.1016/j.ymssp.2025.113239. hal-05229475

HAL Id: hal-05229475

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

Submitted on 29 Aug 2025

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Identification and monitoring of stochastic linear subsystems with unknown local nonlinearities via output injection

Neha Aswalª, Adrien Mélotª, Laurent Mevelª,*, Qinghua Zhangª

ªUniv. Gustave Eiffel, Inria, Cosys-SII, 14S, Campus de Beaulieu, Rennes, France

## Abstract

Most civil and mechanical structures exhibit nonlinear stochastic behaviour, which is difficult to model accurately, but necessary for conventional model-based structural health monitoring techniques. Although various methods have been developed to estimate nonlinear systems, they require information about the external excitation and are susceptible to sensor noise and modelling inaccuracies. This knowledge is challenging to acquire in practice. Hence, this article presents a novel nonlinearity model-agnostic approach to detect damage in mechanical systems with localized nonlinearities. The proposed method utilises output injection to reject the unknown nonlinearities as if they were unknown disturbances. By applying an existing disturbance rejection technique, the need for a priori knowledge about the functional form of nonlinearities is avoided. Besides, a switching strategy is employed to determine the most probable location of the nonlinearities, thereby eliminating the need for a priori knowledge about their location. The method makes use of interacting particle Kalman filter, where the particle filter estimates the parameters (health indices) in order to detect possible damage while the Kalman filter simultaneously estimates the states. The efficiency of the proposed method is demonstrated with the help of numerical experiments on a spring-mass-damper oscillator chain with attached localized nonlinearity. The proposed approach is further validated against experimental data of jointed beams.

**Keywords:** interacting particle Kalman filter; localized nonlinearity; nonlinear vibrations; nonlinearity rejection; structural health monitoring

## 1. Introduction

The importance of regularly monitoring the condition of civil and mechanical structures has been increasingly recognized over the years. To address this need, various structural health monitoring (SHM) algorithms have been developed [1–3], falling into three main categories: data-driven [4, 5], model-driven [6–8], and physics-informed data-driven [9, 10] methods. Data-driven algorithms rely on statistical analysis and machine learning techniques, while model-driven approaches use mathematical and computational models to

<footnote>*Corresponding author; E-mail address: laurent.mevel@inria.fr</footnote>
simulate structural behavior and compare predictions with data. Physics-informed data-driven algorithms have recently emerged in an attempt to bridge the gap between both methodologies by integrating physical principles, e.g. force equilibrium, in machine learning algorithms to enhance predictive accuracy. Although these methodologies seem ideal, they share a significant limitation with purely data-driven methods [3]: the need for an extensive database encompassing numerous potential damage scenarios and environmental conditions to accurately assess structural damage [11]. Acquiring such a comprehensive database can be challenging due to economic or feasibility constraints. Conversely, model-driven SHM methods, despite their computational expense, have been developed to robustly detect structural anomalies, even in the presence of uncertainties.

Model-based SHM usually rely on Bayesian filtering techniques to compare the measured vibrational data (e.g. accelerations, velocities, strain, etc.) observed by sensors to those predicted by the corresponding support numerical model [6, 7]. The Bayesian-filtering based SHM methods require the dynamics of the structure to be represented in a state-space formulation where the unobserved states (e.g. displacements, etc.), propagated in time via the Chapman-Kolmogorov formulation [12], are observed through the known sensor measurement outputs [13]. The equation corresponding to the state propagation in time is called the process model, while the relationship between the states and the measured outputs is defined by the measurement models.

Linear models of system dynamics coupled with Kalman filter [14] have been extensively employed [13], for system identification and for detecting subtle changes indicative of damage. However, linearity is an exception and, in practice, mechanical systems are often composed of an assembly of various subsystems connected by interfaces such as bolts [15], frictional contacts [16, 17], seals [18] or gears [19] which often exhibit nonlinear behaviours even in the absence of damage. Since conventional SHM techniques readily assume linear models [3], they often fall short in accurately monitoring nonlinear systems due to the fact that nonlinear systems exhibit multiple phenomena that have no linear counterparts. For instance, methods like maximum likelihood estimation and expectation maximization often encounter irregular estimation issues caused by nonlinearities. To overcome these challenges, substantial modifications are needed, such as integrating smoothing techniques with particle-based approaches [20].

In this context, various methodologies have been developed in the recent years to monitor the condition and health of structures exhibiting nonlinearities, e.g. methods using nonlinear filters such as the unscented [21], extended [22] Kalman Filters, etc. However, all methods assume a known functional form, or a number of candidate forms of the nonlinearity with parameters to be estimated [23–29]. This induces difficulties, as modelling nonlinear systems is a particularly difficult endeavour, for nonlinear forces often result from a complex interplay between various physical mechanisms and they typically exhibit a high sensitivity to various parameters. Furthermore, methodologies to identify nonlinearities often rely on knowledge about the external excitation [30] or a priori identification of the linear part of the model [31], which is incom-
compatible with operational monitoring. Various methodologies have also been proposed to model nonlinear systems from data. Among them, purely data-based machine learning methods such as Sparse Identification of Nonlinear Dynamics (SINDy) [32] have become popular, although they usually lack physical interpretability and tend to perform poorly outside of their training sets. Cenedese et al. [33, 34] proposed a method to identify reduced order models from measurement data. Although this approach benefits from physical interpretability, it lacks interpretability in the context of SHM because linking a data-driven reduced-order model to damage scenarios in a finite element model of a structure is not straightforward, if at all possible.

These limitations underscore the need for more efficient and robust methods to effectively monitor and diagnose nonlinear structural systems. Addressing this gap, the present study proposes an SHM approach capable of detecting, localizing, and quantifying damage in a nonlinear system, consisting of a linear structure with a local nonlinear element. A key feature of the developed approach is that it does not require any prior knowledge of the nonlinearity, apart from its localized nature. A switching strategy [35–37] is employed to identify the most probable location of the nonlinearity based on likelihoods at each time-step. The need for information on nonlinearity is eliminated by rejecting the nonlinear terms from the state-space formulation through output injection [7, 38]. The reformulated model has been used to estimate the states and structural parameters (health indices) via an interacting particle Kalman filter (IPKF) [7, 39]. It should be noted that in the present methodology, IPKF can be substituted by other filters [40, 41] provided they are associated with the nonlinearity rejection and the unknown position localization mechanisms proposed in this article.

While output injection techniques have traditionally been used to reject inputs considered as perturbations in linear systems, here internal nonlinearities are treated as external perturbations using one such technique, without compromising performance or validity. Furthermore, for simplicity, the external force acting on the structure, also unknown, is assumed as a stationary white Gaussian noise (SWGN) by the algorithm. This assumption can be relaxed to encompass non-stationary excitations, provided that suitable algorithms are used [29, 40–42] to estimate the unknown excitation along with the states.

The paper is structured as follows: the updated state-space equation, after rejection of the unknown nonlinearity, is formulated in Section 2. The switching strategy to find the location of nonlinearity is discussed in Section 3 and Section 4 details the general methodology. Sections 5 and 6 present validations of the approach on numerically and experimentally-generated data set, respectively.

## 2. Nonlinear state-space model formulation - nonlinearity rejection

The present article considers mechanical system with localized nonlinearity, whose dynamics governed by Newton's laws, is described by the following equation,

$$ M\ddot{q}(t) + C\dot{q}(t) + Kq(t) + f_{nl}(q(t), \dot{q}(t)) = f_{ext}(t) \quad (1) $$
where, $\mathbf{q}(t) \in \mathbb{R}^{n_q}$ is the displacement vector of $n_q$ dofs at the time instant $t$, $\dot{\mathbf{q}}(t)$ and $\ddot{\mathbf{q}}(t)$ are the corresponding velocity and acceleration vectors, $\mathbf{M} \in \mathbb{R}^{n_q \times n_q}$ is the mass matrix, and $\mathbf{C} \in \mathbb{R}^{n_q \times n_q}$ is the damping matrix. In addition, $\mathbf{K} \in \mathbb{R}^{n_q \times n_q}$ is the linear stiffness matrix, $\mathbf{f}_{nl}(\mathbf{q}(t), \dot{\mathbf{q}}(t))$ is the vector of nonlinear forces, and $\mathbf{f}_{ext}(t) \in \mathbb{R}^{n_q}$ is the external force acting on the system.

Equation (1) is further written in state-space formulation, as follows,

$$ \dot{\mathbf{x}}(t) = \mathcal{F}\mathbf{x}(t) + \mathcal{B}\mathbf{u}(t) + \mathcal{E}\xi(\mathbf{x}(t)) \quad (2) $$

where, $\mathbf{x} = [\mathbf{q}' \ \dot{\mathbf{q}}']' \in \mathbb{R}^{2n_q}$ is the state vector, with $\bullet'$ indicating the matrix transpose operation, $\mathcal{F}$ is the state matrix, $\mathcal{B}$ is the input matrix, $\mathcal{E}$ maps the effect of nonlinearity on the states. Specifically,

$$ \mathcal{F} = \begin{bmatrix} \mathbf{0}_{n_q \times n_q} & \mathbf{I}_{n_q \times n_q} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix}_{2n_q \times 2n_q}, \quad \mathcal{B} = \begin{bmatrix} \mathbf{0}_{n_q \times n_u} \\ \mathbf{M}^{-1}\tau_b \end{bmatrix}_{2n_q \times n_u} \quad \text{and} \quad \mathcal{E} = \begin{bmatrix} \mathbf{0}_{n_q \times n_f} \\ \mathbf{M}^{-1}\tau_e \end{bmatrix}_{2n_q \times n_f} $$

where $\mathbf{0}$ and $\mathbf{I}$ denote the null and identity matrices, respectively. In addition, $\tau_b\mathbf{u}(t) = \mathbf{f}_{ext}(t)$ and $\tau_e\xi(\mathbf{x}(t)) = \mathbf{f}_{nl}(\mathbf{x}(t))$, where $\tau_b$ and $\tau_e$ are respectively $n_q \times n_u$ and $n_q \times n_f$ boolean selection matrices that define the dofs affected by the external input forces and nonlinear forces. $n_u$ is the number of dofs on which the external force is acting and $n_f$ is the number of non-zero components in the nonlinear force vector $\mathbf{f}_{nl}$, i.e. the dofs that are affected by nonlinearities. In both matrices $\tau_b$ and $\tau_e$, the components are set to one if the corresponding degrees of freedom are affected by external input/nonlinear forces and set to zero otherwise. $\mathbf{u}_{n_u \times 1}$, representing the external input force, is assumed to be SWGN of unknown covariance $(\mathbf{Q}^u)$.

Further, the measurement equation is given by the following,

$$ \mathbf{y}(t) = \mathcal{H}\mathbf{x}(t) + \mathcal{D}\mathbf{u}(t) + \mathcal{J}\xi(\mathbf{x}(t)) + \mathbf{v}(t) \quad (3) $$

where, $\mathbf{y}(t) \in \mathbb{R}^{n_y}$ is the acceleration measurement vector (corresponding to some components of $\ddot{\mathbf{q}}(t)$ measured with conventional accelerometers), $\mathcal{H} \in \mathbb{R}^{n_y \times 2n_q}$ is the output matrix, $\mathcal{D} \in \mathbb{R}^{n_y \times n_u}$ and $\mathcal{J} \in \mathbb{R}^{n_y \times n_f}$ transfer the effect of input and nonlinearity to the measurement outputs, respectively. They are described as, $\mathcal{H} = \tau_y [-\mathbf{M}^{-1}\mathbf{K} \ -\mathbf{M}^{-1}\mathbf{C}]$, $\mathcal{D} = \tau_y [\mathbf{M}^{-1}\tau_b]$, and $\mathcal{J} = \tau_y [\mathbf{M}^{-1}\tau_e]$. $\tau_y$ is an $n_y \times n_q$ Boolean matrix, selecting $n_y$ dofs, describing the available measurement channels. $\mathbf{v} \in \mathbb{R}^{n_y}$ is the measurement uncertainty treated as an SWGN with known covariance of $\mathbf{R}^v$.

Equations (2) and (3) are discretized with a constant sampling period $\Delta t$, such that, at the discrete time instants $k\Delta t$ ($k \in 1, 2, ... N$), the discrete time state $\mathbf{x}_k$ corresponds to $\mathbf{x}(k\Delta t)$. The input $\mathbf{u}(t)$ and $\mathbf{y}(t)$ are discretized accordingly, and the nonlinear term is discretized as

$$ \tau_e \xi_k = \tau_e \xi(\mathbf{x}(k\Delta t)) = \mathbf{f}_{nl}(\mathbf{x}(k\Delta t)) \quad (4) $$

where $\tau_e$ encompasses information regarding the location of nonlinearity, while the discretized nonlinearity is compactly denoted by $\xi_k$. Thus the discretized process and measurement equations are as follows,
$$x_{k+1} = Fx_k + Bu_k + E\xi_k \quad (5)$$

$$y_k = Hx_k + Du_k + J\xi_k + v_k \quad (6)$$

where, $F = e^{\mathcal{F}\Delta t}$, $B = (F - I)\mathcal{F}^{-1}\mathcal{B}$, $E = (F - I)\mathcal{F}^{-1}\mathcal{E}$, $H = \mathcal{H}$, $D = \mathcal{D}$, $J = \mathcal{J}$. It should be noted that $\tau_e$ is housed in $E$ and $J$ (cf. Equations (5) and (6), respectively) and is later dealt with switching strategy described in Section 3. This discretization assumes that the sampling period $\Delta t$ is small enough so that $u(t)$ and $\xi(x(t))$ vary insignificantly during each sampling period.

The state-space equations presented above, whether in continuous time or in discrete time, describe the nominal dynamics behavior of the considered structure. If a damage occurs, the actual dynamics deviates from the nominal behavior. This paper focuses on the detection of damage leading to losses in local stiffness elements. Let $\kappa^{(m)}$ denote the nominal value of the $m^{th}$ stiffness element. At the discrete time instant $k$, the actual stiffness value of this element is equal to $\theta_k^{(m)}\kappa^{(m)}$, with some health index $\theta_k^{(m)} \in [0, 1]$. After a damage implying stiffness loss, $\theta_k^{(m)}$ switches from 1 to another value. The instance of damage occurrence and the health index value after the damage are both unknown to the damage detection algorithm.

By taking into account the evolution of the health indices, the matrices involved in Equations (5) and (6) are time varying. Henceforth, they will be written as $F_k$, $B_k$, $E_k$, and $H_k$ in this paper.

## 2.1. Rejection of nonlinearity

From Equations (5) and (6), it is clear that the process and measurement equations are dependent on the nonlinear term $\xi_k$. This nonlinear term can be rejected out thanks to the redundancy between the process and measurement equations (Equations (5) and (6)). Equation (6) can be re-arranged as: $y_k - H_kx_k - Du_k - J\xi_k - v_k = 0$, which is then inserted into the process equation, yielding

$$ \begin{aligned} x_{k+1} &= F_k x_k + B_k u_k + E_k \xi_k + G_k (y_k - H_k x_k - Du_k - J \xi_k - v_k) \\ &= (F_k - G_k H_k) x_k + (B_k - G_k D) u_k + (E_k - G_k J) \xi_k + G_k y_k - G_k v_k \end{aligned} \quad (7) $$

The matrix $G_k$ is chosen such that $(E_k - G_k J) = 0$, then the nonlinear term ($\xi_k$, as in Equation (4)) is eliminated from the above equation. Thus, $G_k$ is computed as $E_k J^\dagger$, with $\dagger$ denoting matrix Moore-Penrose pseudo-inverse operation $(J^\dagger = (J'J)^{-1}J')$. To ensure feasible and stable state estimation [38], it is assumed that $J$ has a full column rank, implying that the number of measurement outputs is larger than the number of nonlinearity inputs to be rejected. Equation (7) is then re-written as follows

$$x_{k+1} = \tilde{F}_k x_k + \tilde{B}_k u_k + G_k y_k + \tilde{w}_k \quad (8)$$

with more compact notations $\tilde{F}_k = F_k - G_k H_k$, $\tilde{B}_k = B_k - G_k D$, and $\tilde{w}_k = -G_k v_k$.

Although Equation (8) presents $x_k$ in a seemingly linear form, it remains the same nonlinear state as defined in Equation (5). The nonlinearity in Equation (5), which stems from the term $\tau_e \xi_k = \tau_e \xi(x_k) =$
$f_{nl}(x(k\Delta t))$, is masked in Equation (8) by the injected output vector $y_k$. Since $y_k$ is known from sensor measurements, it can be treated as a known input, allowing the application of linear Kalman filter theory.

Moreover, the measurement equation, cf. Equation (6), is decoupled from the nonlinear term $\xi_k$ by projecting the actual measurement ($y_k$) using a transformation matrix $T$, such that $TJ = 0$ [43]. The transformed measurement ($z_k$) equation is then given by the following,

$$z_k = \tilde{\mathbf{H}}_k x_k + \tilde{\mathbf{D}} u_k + \tilde{\mathbf{v}}_k \quad (9)$$

where, $z_k = Ty_k$, $\tilde{\mathbf{H}}_k = TH_k$, $\tilde{\mathbf{D}} = TD$, and $\tilde{\mathbf{v}}_k = Tv_k$.

### 3. Switching strategy - identification of the nonlinearity location

Although output injection avoids the need for a priori assumption about the functional form of the nonlinearity acting on the system, this rejection requires the information regarding the nonlinearity location, which is specified by the Boolean matrix $\tau_e$. To find the location of nonlinearity, it is proposed to try all possible nonlinearity locations and identify the most likely location based on the likelihoods of these locations evaluated together with state and parameter estimation. The likelihood is a conditional probability, whose formulation depends on the Kalman innovation, whose computation is conditioned on the different possible nonlinearity locations. Depending on the likelihood values, the location of nonlinearity is determined as the one associated to the highest likelihood value. Thus, the state and parameter estimates, to be transferred for next iteration, are chosen corresponding to the identified location. This proposition is the essence of switching strategy which has been parallelized in an effort to reduce computational cost incurred from running multiple cases.

With the switching strategy, multiple possible locations are considered by introducing $1 \le h \le n_\tau$ different $\tau_e^h$ models, and the state and the health index estimations corresponding to different models are made in parallel. The most plausible nonlinearity location model is selected on the switch value ($h^* \in \{1, 2, \dots, n_\tau\}$) determined by a likelihood factor described later in Section 4, without assuming any prior probability about the nonlinearity location. Then, the state and parameter estimations associated with the selected nonlinearity location, corresponding to $\tau_e^{h^*}$, are carried on to the next time-step.

We emphasize here that while the methodology can theoretically handle any combination of degrees of freedom – provided that the condition of localized nonlinearity and the disturbance rejection criteria [38] holds – its performance is also dependent on the choice and the number of potential locations considered. This can lead to a prohibitive number of combinations and computational burden in the case of high-dimensional finite-element models. Therefore, any available practical knowledge about likely locations should be used to reduce computational costs, either by assigning priors for the switches, or reducing the number of possible locations. In practical engineering applications, while the exact location of the nonlinearity may be unknown, it is often possible to make informed a priori guesses about several likely locations, e.g. around connections
between mechanical components. This prior knowledge can significantly reduce the number of combinations the algorithm needs to consider. For example, if one assumes that the nonlinear force only acts on a single dof of the numerical model, $\tau_e$ reduces to a vector of size $n_q \times 1$ with $n_f = 1$ and the $n_\tau$ possible models corresponding to the different nonlinearity locations are:

$$[\tau_e^1 \quad \tau_e^2 \quad \dots \quad \tau_e^{n_\tau}] = \begin{bmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & & 1 \end{bmatrix}_{n_q \times n_\tau} \quad (10)$$

The flow of the switching strategy has been described by the pseudo-code, Algorithm 1 given in Appendix A.

## 4. Methodology

The proposed damage detection method for nonlinear systems is based on the estimation of the health indices that are associated with the stiffness elements (described in Section 2) with simultaneous estimation of nonlinearity location. To accomplish this task, following difficulties have been addressed.

1. The considered structure involves unknown nonlinearities. They are rejected via output injection.

2. The locations of the nonlinearities are unknown. To address this difficulty, the likelihoods of the presumed location configurations are computed, and the configuration corresponding to the maximum likelihood value is selected.

3. The estimation of the health indices, involved in the state-space model matrices $F_k$, $B_k$, etc., depends on the unknown state vector $x_k(\theta_k)$. It is addressed by applying IPKF, where the states are estimated using KF which is nested inside the PF that estimates the health indices.

As the actual location of the nonlinearity is unknown, multiple instances of the damage detection algorithm (detailed in the following subsections) are run in parallel, each assuming a particular nonlinearity location. When the correct location is guessed, the state-parameter estimation yields the highest likelihood. The $h^*$ estimated nonlinearity location (NL) model is specified with a Boolean vector $\tau_e^{h^*}$, as one of the columns of the matrix in Equation (10). Here, $h^*$ is the switch value which is to be estimated to compute the NL model at a particular time. At $k^{th}$ time-step, $h$-superscript is assigned to the possible NL models denoted by $\tau_e^h$, until the unknown NL model is estimated as $\tau_{e_k}^{h^*}$. Within a time-step, $h$ has all switch values from 1 to $n_\tau$, since the multiple possible models ($\tau_e^{1:n_\tau}$) are run parallel to one another. For each possible $h^{th}$ model, the IPKF updated with the proposed nonlinearity state-space formulation (named sr-IPKF) is run to estimate the corresponding possible states and parameters ($\theta$). The 'sr' (nonlinearity rejection with
207

switching) in sr-IPKF represents the transformed state-space model and is used to distinguish it from the
existing IPKF method.

208

209

Further, a likelihood corresponding to each $h^{th}$ NL model is computed, which is dependent on the
measurement error and its covariance. From the likelihoods, normalised weights, for each ($h^{th}$) NL model
transition, are determined and the most probable $h^{th}$ model for $k^{th}$ step is selected ($= \tau_{ek}^{h*}$) based on the
probable model transition with maximum weight. Following the determination of the switch value $h_k^*$,
the corresponding states, covariance, and parameters are selected from the possible choices (obtained from
sr-IPKF) and carried onto the next time-step.

210

211

212

213

214

215

The computational time to execute sr-IPKF for multiple NL models is compensated by implementing
parallel execution on multiple workers. In the following, sr-IPKF is discussed with the switching strategy
implementation. Furthermore, the selection criterion to determine the most probable model for the particular
time-step is also described. A workflow showing the general principles of the proposed methodology is
depicted in Figure 1.

216

217

218

219

220

## 4.1. Interacting particle and Kalman filter with nonlinearity rejection - sr-IPKF

221

The proposed methodology estimates the health indices corresponding to all $n_m$ stiffness elements, which are encompassed within the vector $\theta$, wherein each component's value ranges from 0 to 1. Here, 0 and 1 correspond to a 100% damage in the element and perfectly healthy state of the element, respectively. Since $\theta$ and $\mathbf{x}$ estimations are inter-dependent, as state-space matrices require $\theta$, which in turn needs $\mathbf{x}$ estimation, the interacting particle Kalman filter (IPKF) is applied to jointly estimate the health indices ($\theta$) and the state vector ($\mathbf{x}$). Though nonlinear systems are considered in this work, the nonlinearity term is rejected by means of output rejection, resulting in the apparently linear state-space model expressed in Equations (8) and (9), hence KF has been readily applied for state estimation in the present study.

222

223

224

225

226

227

228

229

At the initialization step, for each $h(= 1 : n_{\tau})$ possible NL models, all the health indices are initiated
at 1 (healthy damage-free state) following a normal distribution that is applied to obtain corresponding $N_p$
particle vectors. These $N_p$ vectors, denoted by $\theta^1, \theta^2, ..., \theta^{N_p}$, are arranged in an $n_m \times N_p$ matrix, $\mathbf{\Theta}_0$.
Accordingly $N_p$ different versions of the $\theta$-dependent systems matrices $\tilde{\mathbf{F}}_0, \tilde{\mathbf{H}}_0$, etc., are computed and a
corresponding KF is implemented for state estimation. Therefore, $N_p$ KFs are implemented.

230

231

232

233

234

At the $k^{th}$ time-step, the state estimates and the particles are both updated. Each $j^{th}$ particle is evolved
through a random perturbation around its current position. The evolution of the particles to the current $k^{th}$
time-step is undertaken through Gaussian blurring, $\theta_k^{jh} = \alpha \theta_{k-1}^{jh} + N(\delta\theta_{k-1}, \sigma_{k-1}^\theta)$ [7], with $1 \le j \le N_p$. $\alpha$ is a hyper-parameter, tuning the oscillations in $\theta$ estimation. The propagated particles are passed through
KF for state evolution, where they have been employed to define the system matrices $\tilde{\mathbf{F}}_k^{jh}, \tilde{\mathbf{B}}_k^{jh}, \text{and } \mathbf{E}_k^{jh}$. In
addition, each KF yields a likelihood value, computed from its innovation (also known as prediction error).
The particles evolve based on their estimated likelihood depending on the current output and are further

235

236

237

238

239

240
<figure><img src="image_3.png" /><figcaption>Figure 1: Workflow of the proposed methodology.</figcaption></figure>

re-sampled [44] according to the likelihood values of the corresponding KFs. The whole IPKF algorithm is detailed in what follows.

## Kalman filter - state update

The state ($\mathbf{x}_{k+1|k}^j$) and its covariance ($\mathbf{P}_{k+1|k}^j$) are computed for $h^{th}$-NL models, based on the updated state-space formulations when the inputs $\mathbf{u}_k$ are unknown to the SHM (KF-based) algorithm, given by the following,
$$ \mathbf{x}_{k+1|k}^{j^h} = \tilde{\mathbf{F}}_k^{j^h} \mathbf{x}_{k|k}^{j^h} + \mathbf{G}_k^{j^h} \mathbf{y}_k \quad (11) $$

$$ \mathbf{P}_{k+1|k}^{j^h} = \tilde{\mathbf{F}}_k^{j^h} \mathbf{P}_{k|k}^{j^h} \tilde{\mathbf{F}}_k^{j^h'} + \tilde{\mathbf{B}}_k^{j^h} \mathbf{Q}^u \tilde{\mathbf{B}}_k^{j^h'} + \mathbf{G}_k^{j^h} \mathbf{R}^v \mathbf{G}_k^{j^h'} \quad (12) $$

The actual sensor measurement output, provided to the algorithm, is transformed with $\mathbf{T}^h$ and the innovation is computed as

$$ \epsilon_k^{j^h} = \mathbf{T}^h \mathbf{y}_k - \tilde{\mathbf{H}}_k^{j^h} \mathbf{x}_{k+1|k}^{j^h} \quad (13) $$

Kalman gain is computed as follows,

$$ \begin{aligned} \mathbf{K}_k^{j^h} &= \mathbf{P}_{k+1|k}^{j^h} \tilde{\mathbf{H}}_k^{j^h'} (\tilde{\mathbf{H}}_k^{j^h} \mathbf{P}_{k+1|k}^{j^h} \tilde{\mathbf{H}}_k^{j^h'} + \tilde{\mathbf{D}}^{j^h} \mathbf{Q}^u \tilde{\mathbf{D}}^{j^h'} + \mathbf{T}^h \mathbf{R}^v \mathbf{T}^{h'})^{-1} \\ &= \mathbf{P}_{k+1|k}^{j^h} \tilde{\mathbf{H}}_k^{j^h'} \mathbf{S}_k^{j^h-1} \end{aligned} \quad (14) $$

Based on the Kalman gain, the predicted state and its covariance are updated by the following equations,

$$ \mathbf{x}_{k+1|k+1}^{j^h} = \mathbf{x}_{k+1|k}^{j^h} + \mathbf{K}_k^{j^h} \epsilon_k^{j^h} \quad (15) $$

$$ \mathbf{P}_{k+1|k+1}^{j^h} = (\mathbf{I} - \mathbf{K}_k^{j^h} \tilde{\mathbf{H}}_k^{j^h}) \mathbf{P}_{k+1|k}^{j^h} \quad (16) $$

## Particle filter - particles (health indices) update

The likelihood of each particle is calculated based on the innovation mean and covariance of each KF, assumed to follow normal distribution [7]. $L(\epsilon_k^{j^h} | \theta_k^j)$ is the likelihood for $j^{th}$ particle, which is computed as

$$ L_k^{j^h} = L(\epsilon_k^{j^h} | \theta_k^j) = \frac{1}{\sqrt{2\pi |\mathbf{S}_k^{j^h}|}} e^{-0.5 (\epsilon_k^{j^h})' \mathbf{s}_k^{j^h-1} \epsilon_k^{j^h}} \quad (17) $$

This likelihood information is used to compute the normalized weight for each particle. The weight of $j^{th}$ parameter particle for $h^{th}$ NL model is then given as

$$ w_k^{j^h} = w^h(\theta_k^j) = \frac{L_k^{j^h}}{\sum_{j=1}^{N_p} L_k^{j^h}} \quad (18) $$

The particles, based on their weights, are resampled to replace particles that are outside 90% confidence interval of the $\Theta_k$ using residual resampling method, the details of which are out of scope of this article and can be found in [44, 45]. The likelihoods (cf. Equation (17)) and weights (cf. Equation (18)) computed from this step has been obtained for each of the $n_\tau$ NL models and have been utilized in the following subsection to determine the switch value, in turn estimating the NL model for the $k^{th}$ iteration.
## 4.2. Switch selection

The predicted transformed measurement for each $h^{th}$ NL model, $\bar{z}_k^h$, is obtained by computing the weighted mean of the predicted measurement obtained corresponding to each particle, given as follows,

$$ \bar{z}_k^h = \sum_{j=1}^{N_p} w_k^{j^h} \tilde{\mathbf{H}}_k^{j^h} \mathbf{x}_{k+1|k}^{j^h} \quad (19) $$

Further, the transformed measurement error (innovation) for each $h^{th}$ NL model transition is calculated ($\bar{\varepsilon}_k^h = \mathbf{z}_k - \bar{z}_k^h$). The corresponding innovation covariance is similarly obtained as,

$$ \bar{S}_k^h = \sum_{j=1}^{N_p} w_k^{j^h} \mathbf{S}_k^{j^h} \quad (20) $$

The corresponding likelihood is calculated by the following function [7],

$$ \bar{L}_k^h = \bar{L}(\bar{\varepsilon}_k^h | \tau_{e_k}^h) = \frac{1}{\sqrt{2\pi |\bar{S}_k^h|}} e^{-0.5(\bar{\varepsilon}_k^h)' \bar{S}_k^{h-1} \bar{\varepsilon}_k^h} \quad (21) $$

The probability of $h^{th}$ model being the estimated NL model at $k^{th}$ time-step is computed as the following,

$$ \bar{w}_k^h = \bar{w}(\tau_{e_k}^h) = \frac{\bar{L}_k^h}{\sum_{h=1}^{n_\tau} \bar{L}_k^h} \quad (22) $$

From $\bar{w}_k^h$, the switch value for $k^{th}$ iteration ($h_k^*$) is determined as the index corresponding to the maximum value of the probability of $h^{th}$ model being the estimated NL model at $k^{th}$ time-step.

$$ h_k^* = \operatorname*{argmax}_{1<h<n_\tau} \bar{w}_k^h \quad (23) $$

## 4.3. State and parameter estimation

Following the identification of the NL model at $k^{th}$ time-step as $\tau_{e_k}^{h_k^*}$ from $h_k^*$ determined in Equation (23), the corresponding weights for state and parameter estimation are selected from the overall weights estimated earlier in Section 4.1 (cf. Equation (18)).

$$ w_k^j = w_k^{j^{h_k^*}} \quad (24) $$

Similarly, $\mathbf{x}_{k+1|k+1}^{j} (= \mathbf{x}_{k+1|k+1}^{j^{h_k^*}})$, $\mathbf{\Theta}_k (= \mathbf{\Theta}_k^{h_k^*})$, and $\mathbf{P}_{k+1|k+1}^{j} (= \mathbf{P}_{k+1|k+1}^{j^{h_k^*}})$ are also obtained. These states ($\mathbf{x}_{k+1|k+1}^{j}$), parameters ($\mathbf{\Theta}_k = [\theta^1 \theta^2 \dots \theta^{N_p}]_k$) and covariances ($\mathbf{P}_{k+1|k+1}^{j}$), corresponding to the particles associated to $h_k^*$ model, are carried on to the next time step and the procedure (from Section 4.1) is repeated for further estimations. The particle approximations for the parameter and state means are computed as
$$ \mathbf{x}_{k+1|k+1} = \sum_{j=1}^{N_p} w_k^j \mathbf{x}_{k+1|k+1}^j; \qquad (25) $$

$$ \bar{\theta}_k = \sum_{j=1}^{N_p} w_k^j \theta_k^j \qquad (26) $$

An element $m$ is considered damaged when the corresponding value of the health index $\bar{\theta}_k^{(m)}$ persistently remains below 1, indicating a deviation from the undamaged state. The level at which $\bar{\theta}_k^{(m)}$ value appears to stabilize below 1 provides a quantitative measure of the extent of damage. Due to the fluctuations of the health indices between iterations, distinguishing between damaged and undamaged states can be challenging, potentially leading to false detection with low damage levels. To mitigate this, a threshold should be established; in the examples presented hereafter, a threshold value of 0.8 was deemed satisfactory for accurate damage identification.

The proposed methodology is further described by a pseudo-code given in Algorithm 2 in Appendix A.

## 5. Validation with numerical study

The proposed method simultaneously estimates the location of the nonlinearity as well as the states and parameters of the system. The states are estimated by the KF whereas the parameters ($\theta$) are estimated through PF with the help of 1000 particles ($N_p$). $\alpha$ has been selected as 0.99 for the particle evolution, cf. Section 4.1, for all cases. The selection of $N_p$ and $\alpha$ values has been made empirically, keeping in mind the computational cost and quality of $\theta$ estimates.

### 5.1. Data generation

The numerical experiments are performed on an 8-dof spring-mass-damper system (Figure 2), whose mechanical properties are reported in Table 1. A band-limited white Gaussian noise with root mean square value of 100 N in frequency range [5, 150] Hz is considered to act as external force $f_{ex}(t)$, on all dofs. Linear viscous damping is introduced in the system by considering a stiffness-proportional damping term ($C = \beta K$), with coefficient $\beta = 4.7e^{-3} s^{-1}$ chosen to ensure the damping ratio equals 1% at the frequency of the first linear mode. Although stiffness-proportional damping has been assumed, the proposed methodology is compatible with other damping models as well.

<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic representation of the considered mechanical system with a nonlinear element attached to the eighth dof.</figcaption></figure>
<table><caption>Table 1: Material properties of the mass-spring system.</caption><thead><tr><th>m<sub>1:8</sub> (kg)</th><th>k (N/m)</th><th>k<sub>nl</sub> (N/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>1400</td><td>3e6</td><td>1e9</td></tr></tbody></table>

To validate the methodology, datasets were generated using various functional forms of nonlinearities, with all considered cases detailed in Table 2. Eight models, each assuming the nonlinearity acts on a different degree of freedom, are considered as potential locations within the switching algorithm, as detailed in Equation (10). Furthermore, in this example, damage is modelled as a reduction of stiffness of the linear springs connected to the 5th dof. For each scenario, datasets consisting of the time series of the velocities are obtained by solving the equations of motion of the system using Matlab's ODE45 solver and contaminating them with 1% signal-to-noise ratio (SNR) SWGN to mimic sensor noise.

<table><caption>Table 2: List of considered nonlinear elements.</caption><thead><tr><th>S. no.</th><th>Functional form</th><th>Amplitude (k<sub>nl</sub>/c<sub>nl</sub>)</th><th>Location (dof)</th></tr></thead><tbody><tr><td>1</td><td rowspan="6">k<sub>nl</sub>q<sup>3</sup></td><td rowspan="3">10<sup>12</sup> N/m<sup>3</sup></td><td>1</td></tr><tr><td>2</td><td>3</td></tr><tr><td>3</td><td>6</td></tr><tr><td>4</td><td rowspan="2">10<sup>8</sup> N/m<sup>3</sup></td><td rowspan="3">8</td></tr><tr><td>5</td></tr><tr><td>6</td><td>10<sup>4</sup> N/m<sup>3</sup></td></tr><tr><td>7</td><td>k<sub>nl2</sub>q<sup>2</sup> + k<sub>nl3</sub>q<sup>3</sup></td><td>10<sup>8</sup> N/m<sup>2</sup> - 10<sup>12</sup> N/m<sup>3</sup></td><td></td></tr><tr><td>8</td><td>c<sub>nl</sub>q<sup>2</sup></td><td>10<sup>3</sup> Ns<sup>2</sup>/m<sup>2</sup></td><td></td></tr></tbody></table>

## 5.2. Damage sensitivity

The present methodology aims to determine damage in a nonlinear structure without any knowledge of the nonlinearity and external force (assumed unknown SWGN). Observation output data on all the *dofs* has been provided to sr-IPKF for all numerical studies henceforth. To validate the effectiveness of sr-IPKF, 58 s simulated measurement data at all *dofs* (described in Section 5.1) is provided, where a 60% damage is introduced in the initially undamaged structure at $t = 17$ s. Figure 3 highlights the capability of the proposed methodology in detecting, quantifying and locating damage, as the health index $\theta_5$ associated to the 5th element quickly decreases to a value of $\theta = 0.4$, corresponding to a 60% damage. Meanwhile $\theta$s of the other *dofs* remain at a value of 1. It should be noted that all the health indices ($n_\theta = 8$) are simultaneously estimated, however Figure 3 shows the particle dispersion corresponding to the fifth and the eighth element for better readability of the plot. One can also note that the estimation of $\theta$ over time is stable.

Furthermore, the estimation of the states and the position of the localized nonlinearity is a byproduct of the developed methodology. Figure 4a depicts the comparison of the actual and the estimated velocity at the 8th dof, to which the localised nonlinearity is attached. The estimation is accurate for most parts of the duration, with relative absolute error within a 95% confidence interval (CI) of $(1.47 \times 10^{-4}\%, 8.64\%)$,
<figure><img src="image_5.png" /><figcaption>Figure 3: Particle propagation in time with estimated health indices - sr-IPKF.</figcaption></figure>

cf. Figure 4a. Additionally, as observed in Figure 4b, the sr-IPKF methodology also correctly identifies the location of the nonlinearity as the 8th dof with the highest number of occurrences.

<figure><img src="image_6.png" /><figcaption>Figure 4: Estimation of states (velocity at dof 8) and location of attached nonlinearity.</figcaption></figure>

A comparative study was conducted to demonstrate the efficiency and accuracy for the proposed methodology compared to existing techniques. The estimation of the health index $\theta$ and the state, in this case the velocity of *dof* 8, was chosen as performance metrics. For this study, the measurement (acceleration) data has been simulated for 3072 sample points (≈ 13.4 s) at all *dof* s. Same data length has been used for the following numerical experiment cases. A damage of 60% ($\theta_5$ = 0.4) at element 5. and a nonlinear element
located at the 8th dof.

The three techniques included in the comparison are (1) interacting particle Kalman filter (IPKF), (2) rIPKF (rejection of perturbation with the incorrect location of the nonlinearity assumed at dof 7 instead of dof 8) and (3) sr-IPKF (proposed method with switching strategy). Each of the methods are supplied with the same measurement data. While computationally slower, Figure 5 shows that sr-IPKF is able to determine the correct $\theta$s associated to both damaged and undamaged elements, while both IPKF (1) and rIPKF (2) are not able to estimate the correct health indices. The failure of IPKF is expected due to the lack of nonlinearity rejection. However, method (2) shows that the nonlinearity rejection with an incorrect assumption of nonlinearity location – a common challenge in real-world applications – fails at properly estimating the health indices, as illustrated in Figure 5 (center).

<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of health indices estimation by IPKF, rIPKF (with incorrect assumption for the location of nonlinearity), and proposed method - sr-IPKF; with (d) indicating damage in the element.</figcaption></figure>

In order to further emphasize the efficiency of the methodology, Figure 6 shows the states estimated with all three methods. As expected, the state estimation from IPKF (1) is significantly off from the actual values. Similarly, the state estimation from rIPKF (2) is erroneous due to the incorrect assumption about the location of the nonlinearity. In contrast, the sr-IPKF (3) estimates the states with satisfactory accuracy.

Next, the sensitivity of the proposed method to different levels of damage severity is tested. Four cases have been studied where damage of 20%, 40%, 60% and 80%, respectively have been induced throughout the numerical experiment simulation (≈ 13.4 s constituting 3072 sample points in time for one measurement output). Figure 7 shows that, for all the cases, the health indices are correctly estimated with acceptable accuracy. However, for the 20% damage scenario, the fluctuations of the health indices $\theta$ of the undamaged elements masks the damage estimation. This overlap makes it difficult to distinguish between damaged and undamaged elements. Thus the value of 0.8 for $\theta$ (corresponding 20% damage) has been considered to demarcate the presence of damage in an element.
<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of state estimation (velocity at dof 8) by IPKF, rIPKF (with incorrect assumption for the location of nonlinearity), and proposed method - sr-IPKF.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Sensitivity of the proposed sr-IPKF method to damage severity; with ud and d indicating elements that were undamaged and that incurred damage, respectively.</figcaption></figure>

### 5.3. Effect of various nonlinearities

The proposed method is capable to deal with any type of nonlinearity as it does not assume any functional form. The aim of this section is to demonstrate this by investigating various values of nonlinear coefficients and functional forms with nonlinear stiffness and damping. All cases considered are reported in Table 3. In each case, 60% damage is introduced in the 5th element of the nonlinear structure (shown in Figure 2). The results in terms accuracy of the damage detection are summarized in Table 3. The detection appears accurate, with an accuracy greater than 95% in most cases, with no false positives to report. A slightly smaller accuracy of approximately 90% is observed for cases 6, 7 and 8, corresponding to a relatively mild cubic nonlinearity, a combined quadratic and cubic nonlinearity and a quadratic damping. Although these
results are surprising for case 6 as the nonlinearity is weaker, the accuracy of the health indices $\theta$ estimation is reported as the mean of the last 1000 iterations. Thus, cases 6, 7 and 8 may be able to exhibit better accuracy with more sample points used in the algorithm. The present paper considers 90% accuracy as satisfactory, in order to reliably detect, localize and quantify damage in an element. Thus, only 3072 sample points have been used in this study. Nevertheless, the accuracy ought to be used as a measure to determine the sample points needed to assure smaller fluctuations in the estimation of the health indices $\theta$.

<table><caption>Table 3: Effectiveness of the proposed method against different types of nonlinearity.</caption><thead><tr><th rowspan="2">S. no.</th><th colspan="4">Nonlinearity</th><th colspan="5">Damage</th></tr><tr><th>Functional form</th><th>Amplitude (k<sub>nl</sub>/c<sub>nl</sub>)</th><th>Location (dof)</th><th>Estimation</th><th>Location (element)</th><th>Detection</th><th>Severity</th><th>Accuracy (%)</th><th>False positive</th></tr></thead><tbody><tr><td>1</td><td rowspan="8">k<sub>nl</sub>q<sup>3</sup></td><td rowspan="3">10<sup>12</sup> N/m<sup>3</sup></td><td>1</td><td>✓</td><td rowspan="8">5</td><td>✓</td><td rowspan="8">60%</td><td>94.65</td><td>-</td></tr><tr><td>2</td><td>3</td><td>✓</td><td>✓</td><td>96.13</td><td>-</td></tr><tr><td>3</td><td>6</td><td>✓</td><td>✓</td><td>98.20</td><td>-</td></tr><tr><td>4</td><td rowspan="2">10<sup>8</sup> N/m<sup>3</sup></td><td rowspan="4">8</td><td>✓</td><td>✓</td><td>96.62</td><td>-</td></tr><tr><td>5</td><td>✓</td><td>✓</td><td>95.19</td><td>-</td></tr><tr><td>6</td><td>10<sup>4</sup> N/m<sup>3</sup></td><td>✓</td><td>✓</td><td>92.66</td><td>-</td></tr><tr><td>7</td><td>k<sub>nl2</sub>q<sup>2</sup> + k<sub>nl3</sub>q<sup>3</sup></td><td>10<sup>8</sup> N/m<sup>2</sup> - 10<sup>12</sup> N/m<sup>3</sup></td><td>✓</td><td>✓</td><td>89.09</td><td>-</td></tr><tr><td>8</td><td>c<sub>nl</sub>q<sup>2</sup></td><td>10<sup>3</sup> Ns<sup>2</sup>/m<sup>2</sup></td><td>✓</td><td>✓</td><td>✓</td><td>91.38</td><td>-</td></tr></tbody></table>

## 5.4. Sensor noise sensitivity

Measurement data from real-life sensors are subjected to noise. Hence a sensor noise sensitivity study has been conducted to ensure the robustness of the proposed method to measurement noise. Various SNR levels 1%, 2%, 5% and 10% have been considered in the data obtained by time integration of the equations of motion, with all other parameters kept identical. One can see in Figure 8 that the estimated health indices are close to the actual $\theta$ values and the proposed methodology is mostly insensitive to the different levels of sensor noise.

<figure><img src="image_10.png" /><figcaption>Figure 8: Sensitivity of the proposed sr-IPKF method to measurement noise; with ud and d indicating elements that were undamaged and that incurred damage, respectively.</figcaption></figure>
# 6. Application to a nonlinear jointed system

The acceleration data for the experimental study has been obtained from the dataset provided by [46, 47]. The experiment has been conducted on a modified half Brake-Reuß beam which consists of two 304 stainless steel beams joined together by a 3-bolt lap joint. The overall length of the beam is 0.72 m with a width of 2.54 · 10⁻² m and a depth of 1.27 · 10⁻² m. The details of the experimental setup and procedure to obtain acceleration data are described in Section 2 of Chen et al. [46] and should be referred for more information on the experiment. Although the experimental set-up consists of two bolted beams, in the present study, a linear finite element model of a 2-dimensional (2D) Euler-Bernoulli free-free beam has been chosen as the numerical support model for the modified half Brake-Reuß beam [46, 47]. This simple model is motivated by the small number of available sensors (three accelerometers) since the nonlinear rejection algorithm requires more sensors than the number of nonlinear degrees of freedom to be rejected. Figure 6 shows the location of measurement channels (accelerometers measuring acceleration in translational *dofs* of nodes 2, 7, and 8) and element-wise numerical break-up of the beam’s support model, used for the proposed methodology. The numerical support model has 8 elements and 9 nodes with 2-*dofs* at each node.

<figure><img src="image_11.png" /><figcaption>Figure 9: Numerical representation of the nonlinear system, modified half Brake-Reuß beam [46, 47].</figcaption></figure>

The modulus of elasticity and density were chosen as 195 GPa and 7800 kg/m³, respectively, in accordance with the material properties of the 304 stainless steel used in the experiment. Furthermore, a stiffness-proportional damping model was applied with a damping ratio of 0.1%, as observed in [47] for most acceleration amplitudes. Acceleration data from the shaker ringdown test [46, 47] was employed, where the beam was excited with a shaker until the steady-state response was reached and then allowed to decay freely by removing the external force. Data corresponding to the transient vibration state (free decay) was used for this study, comprising 3072 sample points at a rate of 6400 Hz (iterations 193198 to 196269). It is important to note that the structure is undamaged [46], thus the effectiveness of the proposed method is validated by correctly localizing the nonlinearity (at node 5, as shown in Figure 9) and ensuring the health indices (θ) remain above 0.8 (identified as the threshold value for the presence of damage in Section 5.2).

Similar to numerical studies, 1000 particles are used for the estimation of the health indices with $\alpha = 0.99$. Without loss of generality, the localized nonlinearity is assumed to be acting on any of the translational *dofs* of the numerical support model in order to reduce the computational effort. Hence, there are $n_{\tau} = 9$
possible $\tau_e$ models that are described as follows,

$$ \begin{bmatrix} \tau_e^1 & \tau_e^2 & \dots & \tau_e^9 \end{bmatrix} = \begin{bmatrix} 1 & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \\ 0 & 0 & \dots & 0 \end{bmatrix}_{18 \times 9} \qquad (27) $$

As shown in Figure 10a, the measurement estimation exhibits some deviations from the actual measurements as the algorithm searches for the optimum value of the nonlinearity location, especially at the initial stages. These deviations reduce over time as the computation of location of nonlinearity, $\theta$ and state and innovation covariances improve, reaching relative absolute error within a 95% CI of $(9.7 \cdot 10^{-5}\%, 3.86\%)$ – calculated for acceleration data from 0.2 s (relative time) onwards. Figure 10b shows the spread of $\theta_{1:8}$ estimation, with the minimum, maximum, and mean $\theta$ estimate per iteration emphasized. One can see that temporal mean for $\theta$ ranges from 0.923 in the 8th element to 1.002 in the 1st element, both of which exceed the damage detection threshold of $\theta = 0.8$ (cf. Section 5.2). Additionally, Figure 10c shows that the proposed method is able to successfully identify the location of the nonlinearity, a crucial criterion for validating the methodology.

Since acceleration data with damage in the beam is not available, the algorithm's ability to estimate damage has been tested by monitoring whether $\theta_8$ estimate, initiated at an elevated value (1.3, corresponding to elastic modulus of 253.5 GPa) returns back to its actual value ($\approx 1$). The algorithmic parameters have been kept the same as the previous case, except the acceleration data. The subsequent transient vibration data, corresponding to iterations 196270 to 199341 of the dataset, has been used. The location of nonlinearity, determined previously at the 5th node of the FE model, has also been provided.

Figure 11a shows that proposed method is able to determine the correct value of $\theta_8$ with an accuracy of 97.2% (calculated over the last 1000 iterations). One can also note that, the fluctuation of the other health indices $\theta$ decrease once the initially incorrect $\theta_8$ reaches its correct value, all the while remaining well above the threshold of 0.8. Furthermore, measurement estimation is observed to be fairly accurate (Figure 11b), when compared to the experimental acceleration data, with relative absolute error within a 95% CI of $(1.142 \cdot 10^{-6}\%, 2.73\%)$ – calculated for the entire duration of the acceleration data.

The performance of the proposed methodology could be improved if more measurement channels were available during the experiment, since the output injection technique is limited by the number of measurement channels, which should be greater than the number of rejected inputs, in this case, the *dofs* affected by the nonlinear forces. In this study, although the nonlinearity is localized, the ratio between the num-
<figure><img src="image_12.png" /><figcaption>(a) Measurement estimation</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Health index estimation</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Nonlinearity location estimation</figcaption></figure>

<caption>Figure 10: Estimation of measurement (acceleration at dof 15), health index, and location of nonlinearity.</caption>

ber of nonlinear degrees of freedom and measurement channels is only 1:3. While this satisfies the main assumptions of the methodology, increasing the number of sensors could enhance accuracy. Another way of improving the results, should more measurement channels be available, would be to allow the localized nonlinearity to span several *dof s* in matrix *T*<sub>e</sub>, to better approximate the 3-bolt lap joint, and utilize a more realistic FE model.

## 7. Conclusion

With model-based health monitoring techniques, it is important to have an accurate numerical support model for correct damage detection. However it is cumbersome to know or estimate the exact behaviour of a structure, especially when it exhibits nonlinear behaviours, as is the case for most of the civil and mechanical
<figure><img src="image_15.png" /><figcaption>Figure 11: Estimation of measurement (acceleration at dof 15), and health index; with (d) indicating damage in the element.</figcaption></figure>

structures. The present article proposes a novel methodology which is able to detect, locate and quantify
damage in inherently nonlinear structures with unknown localized nonlinearities.

The proposed method uses output injection technique to update the state-space formulation in order to reject the nonlinear terms from both the process and measurement equations and a switching algorithm to ensure the methodology remains model-agnostic with respect to the nonlinearity. The switching strategy works by introducing multiple possible models, each associated to a possible location of the nonlinearity. At each iteration, the most likely nonlinearity location as well as the system’s states and parameters (health indices) are estimated, based on the data from the previous iteration and measurement data from the present time-step via an interacting particle Kalman filter. The proposed method has been validated against numerically-generated data of an 8-degree-of-freedom spring-mass-damper oscillator chain and experimental data of a modified half Brake–Reuß beam, respectively. The proposed method proved efficient for various types and location of nonlinearities, damage levels and measurement noise levels.

Future work could investigate the scalability of the proposed approach to a larger nonlinear structure and improved performances with smaller rejected localised nonlinearity-to-measurement ratio for wider application to various nonlinear systems. In addition, the current methodology explores each possible nonlinearity position individually, which may act as a limitation due to increased computational efforts. Future improvements could aim at reducing this combinatorial cost by embedding more flexible or adaptive model selection strategies. While this study focuses on damage detection in linear subsystems by rejecting nonlinearities, extending the approach to identify and characterize these nonlinear effects represents another promising research direction which could bring improvement in structural health monitoring and its applications.
456

Acknowledgement

457

This project has received funding in part from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 899546.

458

459

References

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

[1] Hoon Sohn, Charles R Farrar, Francois M Hemez, Devin D Shunk, Daniel W Stinemates, Brett R Nadler, and Jerry J Czarnecki. A review of structural health monitoring literature: 1996–2001. *Los Alamos National Laboratory, USA*, 1:16, 2003.

[2] Deepam Goyal and BS Pabla. The vibration monitoring methods and signal processing techniques for structural health monitoring: a review. *Archives of Computational Methods in Engineering*, 23:585–594, 2016.

[3] Rongrong Hou and Yong Xia. Review on the new development of vibration-based damage identification for civil engineering structures: 2010–2019. *Journal of Sound and Vibration*, 491:115741, 2021.

[4] Charles R Farrar and Keith Worden. *Structural health monitoring: a machine learning perspective*. John Wiley & Sons, 2012.

[5] Mohsen Azimi, Armin Dadras Eslamlou, and Gokhan Pekcan. Data-driven structural health monitoring and damage detection through deep learning: State-of-the-art review. *Sensors*, 20(10):2778, 2020.

[6] Antonios Kamariotis, Eleni Chatzi, and Daniel Straub. Value of information from vibration-based structural health monitoring extracted via bayesian model updating. *Mechanical Systems and Signal Processing*, 166:108465, 2022.

[7] Subhamoy Sen, Neha Aswal, Qinghua Zhang, and Laurent Mevel. Structural health monitoring with non-linear sensor measurements robust to unknown non-stationary input forcing. *Mechanical Systems and Signal Processing*, 152:107472, 2021.

[8] Suzana Ereiz, Ivan Duvnjak, and Javier Fernando Jiménez-Alonso. Review of finite element model updating methods for structural applications. In *Structures*, volume 41, pages 684–723. Elsevier, 2022.

[9] Zhiming Zhang and Chao Sun. Structural damage identification via physics-guided machine learning: a methodology integrating pattern recognition with finite element model updating. *Structural Health Monitoring*, 20(4):1675–1688, 2021.

[10] TG Ritto and FA Rochinha. Digital twin, physics-based model, and machine learning applied to damage detection in structures. *Mechanical Systems and Signal Processing*, 155:107614, 2021.

[11] Onur Avci, Osama Abdeljaber, Serkan Kiranyaz, Mohammed Hussein, Moncef Gabbouj, and Daniel J Inman. A review of vibration-based damage detection in civil structures: From traditional methods to machine learning and deep learning applications. *Mechanical systems and signal processing*, 147:107077, 2021.

[12] John Frank Charles Kingman. On the chapman—kolmogorov equation. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 276(1259):341–369, 1974.

[13] Zhe Chen et al. Bayesian filtering: From kalman filters to particle filters, and beyond. *Statistics*, 182(1):1–69, 2003.

[14] Simon Haykin. *Kalman filtering and neural networks*. John Wiley & Sons, 2004.

[15] H. Festjens, G. Chevallier, and J.L. Dion. Nonlinear model order reduction of jointed structures for dynamic analysis. *Journal of Sound and Vibration*, 333(7):2100–2113, 2014. ISSN 0022-460X. doi: https://doi.org/10.1016/j.jsv.2013.11.039. URL https://www.sciencedirect.com/science/article/pii/S0022460X13009930.

[16] M. Krack, L. Salles, and F. Thouverez. Vibration prediction of bladed disks coupled by friction joints. *Archives of Computational Methods in Engineering*, 24(3):589–636, 2017. ISSN 1886-1784. doi: 10.1007/s11831-016-9183-2. URL https://doi.org/10.1007/s11831-016-9183-2.

[17] Mertol Tüfekci, John P. Dear, and Loïc Salles. Forced vibration analysis of beams with frictional clamps. *Applied Mathematical Modelling*, 128:450–469, 2024. ISSN 0307-904X. doi: https://doi.org/10.1016/j.apm.2024.01.031. URL https://www.sciencedirect.com/science/article/pii/S0307904X24000313.

[18] Mertol Tüfekci, Fadi El Haddad, Loïc Salles, Richard Setchfield, and Ludovic Renson. Effects of the Seal Wire on the Nonlinear Dynamics of the Aircraft Engine Turbine Blades. *Journal of Engineering for Gas Turbines and Power*, 145 (12):121004, 10 2023. ISSN 0742-4795. doi: 10.1115/1.4063413. URL https://doi.org/10.1115/1.4063413.

[19] Adrien Mélot, Joël Perret-Liaudet, and Emmanuel Rigaud. Vibro-impact dynamics of large-scale geared systems. *Nonlinear Dynamics*, December 2022. ISSN 1573-269X. doi: https://doi.org/10.1007/s11071-022-08144-5.

[20] Thomas B Schön, Adrian Wills, and Brett Ninness. System identification of nonlinear state-space models. *Automatica*, 47(1):39–49, 2011.

[21] Rodrigo Astroza, Hamed Ebrahimian, Yong Li, and Joel P Conte. Bayesian nonlinear structural fe model and seismic input identification for damage assessment of civil structures. *Mechanical Systems and Signal Processing*, 93:661–687, 2017.

[22] Jann N. Yang, Silian Lin, Hongwei Huang, and Li Zhou. An adaptive extended kalman filter for structural damage identification. *Structural Control and Health Monitoring*, 13(4):849–867, 2005. doi: 10.1002/stc.84.

[23] Hamed Ebrahimian, Rodrigo Astroza, Joel P Conte, and Raymond A de Callafon. Nonlinear finite element model updating for damage identification of civil structures using batch bayesian estimation. *Mechanical Systems and Signal Processing*, 84:194–222, 2017.

[24] Hamed Ebrahimian, Rodrigo Astroza, Joel P Conte, and Costas Papadimitriou. Bayesian optimal estimation for output-only nonlinear system and damage identification of civil structures. *Structural Control and Health Monitoring*, 25(4): e2128, 2018.
[25] Sha Wei, Xu-Long Li, Hu Ding, and Li-Qun Chen. A two-step method to locate multiple local nonlinearities. *Mechanical Systems and Signal Processing*, 202:110712, 2023.

[26] Patrick Schrangl, Pavlo Tkachenko, and Luigi Del Re. Iterative model identification of nonlinear systems of unknown structure: Systematic data-based modeling utilizing design of experiments. *IEEE Control Systems Magazine*, 40(3):26-48, 2020.

[27] Chao Huang. Identification of polynomial nonlinear systems based on center manifold. *Automatica*, 159:111349, 2024.

[28] Xiang Liu and Guoping Cai. A novel time-domain approach for identifying nonlinear structural dynamical system with explicit model based on observer/kalman filter identification method. *Journal of Sound and Vibration*, page 118363, 2024.

[29] Tianhao Yu, Zhiheng Wang, and Jingfeng Wang. An iterative augmented unscented kalman filter for simultaneous state-parameter-input estimation for systems with/without direct feedthrough. *Mechanical Systems and Signal Processing*, 205:110793, 2023.

[30] Jia He, Xiaoxiong Zhang, Mengchen Qi, and Bin Xu. Model-free identification of nonlinear restoring force with modified observation equation. *Applied Sciences*, 9(2):306, 2019.

[31] Ying Lei, Xiongjun Yang, Jinshan Huang, Fubo Zhang, and Lijun Liu. Identification of model-free hysteretic forces of magnetorheological dampers embedded in buildings under unknown excitations using incomplete structural responses. *Structural Control and Health Monitoring*, February 2021. doi: 10.1002/stc.2715.

[32] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. *Proceedings of the National Academy of Sciences*, 113(15):3932-3937, 2016. doi: 10.1073/pnas.1517384113. URL https://www.pnas.org/doi/abs/10.1073/pnas.1517384113.

[33] M. Cenedese, J. Axàs, H. Yang, M. Eriten, and G. Haller. Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 380(2229):20210194, June 2022. doi: 10.1098/rsta.2021.0194. URL https://doi.org/10.1098/rsta.2021.0194. Publisher: Royal Society.

[34] Mattia Cenedese, Joar Axås, Bastian Bäuerlein, Kerstin Avila, and George Haller. Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds. *Nature communications*, 13(1):872, 2022. URL https://doi.org/10.1038/s41467-022-28518-y.

[35] KP Murphy. Switching kalman filters. *Technical Reports*, 1998.

[36] Yi-Ming Zhang, Hao Wang, Yu Bai, Jian-Xiao Mao, Xiang-Yu Chang, and Li-Bin Wang. Switching bayesian dynamic linear model for condition assessment of bridge expansion joints using structural health monitoring data. *Mechanical Systems and Signal Processing*, 160:107879, 2021.

[37] Neha Aswal, Subhamoy Sen, and Laurent Mevel. Switching kalman filter for damage estimation in the presence of sensor faults. *Mechanical Systems and Signal Processing*, 175:109116, 2022.

[38] Qinghua Zhang and Liangquan Zhang. State estimation for stochastic time varying systems with disturbance rejection. *IFAC-PapersOnLine*, 51(15):55-59, 2018.

[39] Rickard Karlsson, Thomas Schon, and Fredrik Gustafsson. Complexity analysis of the marginalized particle filter. *IEEE Transactions on Signal Processing*, 53(11):4408-4411, 2005.

[40] Ying Lei, Dandan Xia, Kalil Erazo, and Satish Nagarajaiah. A novel unscented kalman filter for recursive state-input-system identification of nonlinear systems. *Mechanical Systems and Signal Processing*, 127:120-135, 2019.

[41] Ying Lei, Yongqiang Jiang, and Zhiqian Xu. Structural damage detection with limited input and output measurement signals. *Mechanical systems and signal processing*, 28:229-243, 2012.

[42] Subhamoy Sen, Antoine Crinière, Laurent Mevel, Frédéric Cérou, and Jean Dumoulin. Seismic-induced damage detection through parallel force and parameter estimation using an improved interacting particle-kalman filter. *Mechanical Systems and Signal Processing*, 110:231-247, 2018.

[43] Eshwar Kuncham, Neha Aswal, Subhamoy Sen, and Laurent Mevel. Bayesian monitoring of substructures under unknown interface assumption. *Mechanical Systems and Signal Processing*, 193:110269, 2023.

[44] Jun S Liu and Rong Chen. Sequential monte carlo methods for dynamic systems. *Journal of the American statistical association*, 93(443):1032-1044, 1998.

[45] Tiancheng Li, Miodrag Bolic, and Petar M Djuric. Resampling methods for particle filtering: classification, implementation, and strategies. *IEEE Signal processing magazine*, 32(3):70-86, 2015.

[46] Wei Chen, Debasish Jana, Aryan Singh, Mengshi Jin, Mattia Cenedese, Giancarlo Kosova, Matthew RW Brake, Christoph W Schwingshackl, Satish Nagarajaiah, Keegan J Moore, et al. Measurement and identification of the nonlinear dynamics of a jointed structure using full-field data, part i: Measurement of nonlinear dynamics. *Mechanical Systems and Signal Processing*, 166:108401, 2022.

[47] Mengshi Jin, Giancarlo Kosova, Mattia Cenedese, Wei Chen, Aryan Singh, Debasish Jana, Matthew RW Brake, Christoph W Schwingshackl, Satish Nagarajaiah, Keegan J Moore, et al. Measurement and identification of the nonlinear dynamics of a jointed structure using full-field data; part ii-nonlinear system identification. *Mechanical Systems and Signal Processing*, 166:108402, 2022.

## Appendix A. Pseudo-codes for proposed algorithm

The switching strategy utilized in the present manuscript has been represented by the pseudo-code given in Algorithm 1.
<caption>**Algorithm 1** Pseudo-code for switching strategy.</caption>

1: **procedure** SWITCHING STRATEGY(**y**<sub>k</sub>, **Q**<sup>u</sup>, **R**<sup>v</sup>) ▷ Input measurement and process and measurement noise covariances
2: At k = 0, initialize {**x**<sup>h</sup><sub>0|0</sub>} and {**P**<sup>h</sup><sub>0|0</sub>} ▷ Initialization
3: **for** k<sup>th</sup> time step **do**
4: **for** all h = 1 to n<sub>τ</sub> **do**
5: Select τ<sub>e</sub><sup>h</sup> corresponding to h<sup>th</sup> model
6: Compute **J**<sup>h</sup>, **T**<sup>h</sup> corresponding to τ<sub>e</sub><sup>h</sup>
7: Transform measurement as **T**<sup>h</sup>**y**<sub>k</sub>
8: **procedure** IPKF(**T**<sup>h</sup>**y**<sub>k</sub>, **Q**<sup>u</sup>, **R**<sup>v</sup>) ▷ Input transformed measurement and process and measurement noise covariances
9: **end procedure** Output: predicted transformed measurement (H̃<sub>k</sub><sup>j</sup> x<sub>k+1|k</sub><sup>j</sup>), S<sub>k</sub><sup>j</sup>, and w<sub>k</sub><sup>j</sup>
10: Compute weighted mean of the predicted transformed measurement obtained corresponding to each particle, z̄<sub>k</sub><sup>j</sup>
11: Compute transformed measurement error (innovation) for each h<sup>th</sup> NL model transition, ε̃<sub>k</sub><sup>j</sup>
12: Compute corresponding innovation, S̄<sub>k</sub><sup>j</sup> and likelihood L̄<sub>k</sub><sup>j</sup>
13: **end for**
14: Estimate probability of h<sup>th</sup> model being the nonlinear model, w̄<sub>k</sub><sup>j</sup>
15: Switch value for k<sup>th</sup> iteration, h<sub>k</sub><sup>*</sup> (= argmax<sub>1&lt;h&lt;n<sub>τ</sub></sub> w̄<sub>k</sub><sup>j</sup>)
16: Update state and parameters corresponding to h<sub>k</sub><sup>*</sup>
17: **end for**
18: **end procedure**

The proposed methodology has been described by the following pseudo-code, cf. Algorithm 2.
<caption>**Algorithm 2** Pseudo-code for the proposed methodology - srIPKF.</caption>

1: **procedure** SR-IPKF(y<sub>k</sub>, Q<sup>u</sup>, R<sup>v</sup>) ▷ Input measurement and process and measurement noise covariances
2: At k = 0, initialize state {x<sup>j</sup><sub>0|0</sub>}, state covariance {P<sup>j</sup><sub>0|0</sub>}, parameters {θ<sub>0</sub>} and particles, {Θ<sub>0</sub>} ▷ Initialization
3: **for** k<sup>th</sup> time step **do**
4: Evolve particles to current time step, θ<sup>j</sup><sub>k-1</sub><sup>h</sup> → θ<sup>j</sup><sub>k</sub>
5: **procedure** SWITCHING STRATEGY(y<sub>k</sub>, Q<sup>u</sup>, R<sup>v</sup>) ▷ Parallel loop
6: **for** all h = 1 to n<sub>τ</sub> **do**
7: Select τ<sup>h</sup><sub>e</sub> corresponding to h
8: Compute J<sup>h</sup>, T<sup>h</sup> corresponding to τ<sup>h</sup><sub>e</sub>
9: Transform measurement as T<sup>h</sup>y<sub>k</sub>
10: **procedure** INTERACTING PARTICLE-KALMAN FILTER(T<sup>h</sup>y<sub>k</sub>, Q<sup>u</sup>, R<sup>v</sup>) ▷ See Section 2.1
11: **for** <each particle θ<sup>j</sup><sub>k</sub> ∈ Θ<sup>j</sup><sub>k</sub> > **do**
12: Define system matrices F<sup>j</sup><sub>k</sub>, B̃<sup>j</sup><sub>k</sub>, E<sup>j</sup><sub>k</sub>, H̃<sup>j</sup><sub>k</sub> and D̃<sup>j</sup>
13: **procedure** KALMAN FILTER ▷ See Equations (11) and (12)
14: Predict state (x<sup>j</sup><sub>k+1|k</sub>) and its covariance (P<sup>j</sup><sub>k+1|k</sub>) ▷ See Equations (11) and (12)
15: Compute innovation (ε<sup>j</sup><sub>k</sub>) and Kalman gain (K<sup>j</sup><sub>k</sub>) ▷ See Equations (13) and (14)
16: Update predicted state and its covariance: x<sup>j</sup><sub>k+1|k+1</sub> and P<sup>j</sup><sub>k+1|k+1</sub> ▷ See Equations (15) and (16)
17: **end procedure**
18: Compute likelihood of j<sup>th</sup> particle, L<sup>j</sup><sub>k</sub> ▷ See Equation (17)
19: **end for**
20: Compute the weight of j<sup>th</sup> parameter particle for h<sup>th</sup> NL model, w<sup>j</sup><sub>k</sub> ▷ See Equation (18)
21: **end procedure** Output: predicted transformed measurement (H̃<sup>j</sup><sub>k</sub> x<sup>j</sup><sub>k+1|k</sub>), S<sup>j</sup><sub>k</sub>, and w<sup>j</sup><sub>k</sub>
22: Compute z̃<sup>j</sup><sub>k</sub> and ε̃<sup>j</sup><sub>k</sub> ▷ See Equations (19)
23: Compute corresponding innovation, S̃<sup>j</sup><sub>k</sub> and likelihood L̃<sup>j</sup><sub>k</sub> ▷ See Equations (20) and (21)
24: **end for**
25: Estimate probability of h<sup>th</sup> model being the nonlinear model, w̃<sup>j</sup><sub>k</sub> ▷ See Equation (22)
26: Estimate Switch value for k<sup>th</sup> iteration, h<sup>*</sup><sub>k</sub> (= argmax<sub>1&lt;h&lt;n<sub>τ</sub></sub> w̃<sup>j</sup><sub>k</sub>) ▷ See Equation (23)
27: **end procedure**
28: Following the identification of the NL model at k<sup>th</sup> time-step as τ<sup>h</sup><sub>e</sub><sup>*</sup>, ▷ See section 4.3
29: Weights (w<sup>j</sup><sub>k</sub>), state (x<sup>j</sup><sub>k+1|k+1</sub>), state covariance (P<sup>j</sup><sub>k+1|k+1</sub>), and particles (Θ<sub>k</sub>) are selected corresponding to h<sup>*</sup><sub>k</sub>
30: Compute final state (x<sub>k+1|k+1</sub>) and parameter (θ̃<sub>k</sub>) estimates for k<sup>th</sup> time-step ▷ See Equations (25) and (26)
31: **end for**
32: **end procedure**