A Discontinuous Unscented Kalman Filter for Non-Smooth Dynamic Problems

1. Introduction

The increasing availability of dense and heterogeneous sensor information has allowed for condition assessment and robust diagnostics of linear and non-linear engineered systems across diverse domains, including the civil, mechanical, and aerospace fields (Kumar and Crassidis, 2007; Worden et al., 2008; Farrar and Worden, 2012). Of particular interest for a number of specialized implementations, including that of robust diagnostics and control, are systems that extend beyond the linear range, commonly attained in response to extreme or unusual loads. Such loads may induce behavior that is non-linear and potentially non-smooth, as in the case of plasticity (Smyth et al., 1999; Ebrahimian et al., 2017), impact (Wriggers, 1991), fracture (Kakouris and Triantafyllou, 2017), and sliding (Giannakopoulos, 1989). The adequate modeling of such systems may be achieved by application of offline (Papadimitriou and Papadioti, 2013; Au and Zhang, 2016; Zhang and Au, 2016) or of (near) “real-time” identification schemes, the latter typically relying on adoption of Bayesian-type filters (Kalman, 1963; Ljung and Glad, 1994). The task of online identification is a non-trivial one, especially when conjoined with identification of system parameters, which are often not known a priori or are highly uncertain (Astroza et al., 2017).

When both unmeasured system states and system parameters are to be estimated, a so-called problem of joint state and parameter identification is posed, often expressed as a non-linear system identification case (Chatzi and Smyth, 2014). In previous years, online non-linear estimation was for the most part achieved by means of the Extended Kalman Filter (EKF) (Mariani and Corigliano, 2005; Ding et al., 2014; Ebrahimian et al., 2015). However, more recently, methods which avoid linearization of the state equations, such as the Ensemble Kalman Filter (EnKF) (Huang et al., 2017), the Unscented Kalman Filter (UKF) (Julier and Uhlmann, 1997; Chatzi et al., 2010; Omrani et al., 2013; Al-Hussein and Haldar, 2015), and the Particle Filter (Chatzi and Smyth, 2012; Eftekhar Azam et al., 2012), have gained in popularity due to their flexibility in treating non-linear dynamics. Of the aforementioned techniques, the UKF in particular employs a reduced number of particles, termed the Sigma-Points, maintaining a rapid and online estimation. This is the main driver behind its selection as the method of choice in this paper.

Kalman-type, and Bayesian filters in general, place a fundamental assumption on the dynamic states and the time-invariant parameters of a system being fully observable (Kalman, 1963; Hermann and Krener, 1977; Diop and Fliess, 1991) and identifiable (Walter, 1982). While this holds true for smooth systems, the same does not apply for their non-smooth counterpart, which pertains to systems that are described by non-differentiable state-space equations (Chatzis et al., 2014; Olivier and Smyth, 2017a). Nonetheless, the simulation and tracking of non-smooth systems is essential for numerous engineering problems, since these are by default tied to manifestations of damage and failure.

Previous work of the authors (Chatzis et al., 2014) overviews the classification of systems in accordance with their observability and identifiability properties. Violation of this property forms a salient obstacle for Bayesian-type online identification algorithms, which are expected to diverge for unobservable states or parameters (Liu et al., 1996), which is naturally amplified for non-linear systems. The methods presented in Chatzis et al. (2014) may be implemented to infer the observability and identifiability of a system’s states and parameters, and the way in which these evolve during a system’s response under specific loads. This information may then be seeded into the online estimators (filters) in real-time for ensuring convergence, thereby improving estimation.

In this work, and following original developments introduced in Chatzis et al. (2017), a modified version of the standard UKF is proposed. The key to the formulation lies in adoption of a modular state-space formulation, evaluation of the observability within each time step of the analysis, and selection of an appropriate sub-space of the full state vector to be used by the UKF. The method is termed the Discontinuous D– modification to the UKF, i.e., the DUKF. Examples are presented illustrating the performance of the method for models used in the context of plasticity, which are however general enough to be applied in several other applications of non-smooth problems. The examples reveal a consistently superior performance of the D– modification further highlighting the effects of the special observability properties of non-smooth problems. The proposed alternative opens up the way for robust online tracking and control of a variety of engineered systems including rocking (Chatzis and Smyth, 2012a,b, 2013; Greenbaum et al., 2015), energy (Alavi et al., 2015), and biological systems (Villaverde et al., 2016; Villaverde and Banga, 2017).

2. Non-Smooth Dynamical Systems

A non-linear system with state variables x_t, time-invariant parameters θ, known input vector u, and measurement vector y can in general be described by the following system of equations:

${\dot{\mathbf{\text{x}}}}_t=E\left({\mathbf{\text{x}}}_t,\mathrm{\theta },\mathbf{\text{u}}\right),\dot{\mathrm{\theta }}=0,\mathbf{\text{y}}=G\left({\mathbf{\text{x}}}_t,\mathrm{\theta },\mathbf{\text{u}}\right)$ (1)

where E and G designate the non-linear state-space and measurement functions, respectively. For uncertain systems, i.e., systems whose time-invariant parameters are uncertain or unknown, the above problem may be recast into one of joint state and parameter identification. In this case, the state-space and measurement equations of formulation (1) may be written in an augmented form by introducing the state vector x = [x_t, θ]:

$\dot{\mathbf{\text{x}}}=e\left(\mathbf{\text{x}},\mathbf{\text{u}}\right),\mathbf{\text{y}}=g\left(\mathbf{\text{x}},\mathbf{\text{u}}\right)$ (2)

In the latter representation, one treats both the dynamic states and the parameters of the system as states of the augmented system. A dynamical system is further characterized as analytic or smooth, when the state-space equation (2) are continuous and infinitely differentiable. Very often, however, the state-space equations of physical models may not be analytic, either due to discontinuities in the state-space equation or in their derivatives. It should also be highlighted that smoothness requires that the equations are infinitely differentiable through at least all the realizations of the states encountered during the trajectory of the system. In this paper, we deal with models for which the state-space equations are continuous, but not differentiable, and whose state-space equations can be separated into smooth, i.e., continuous and infinitely differentiable, branches of the form:

$\dot{\mathbf{\text{x}}}=\left\{\begin{array}{c}\hfill e_1\left(\mathbf{\text{x}}\right),\mathbf{\text{x}}\in R_1^n\hfill \\ \hfill ⋮\hfill \\ \hfill e_i\left(\mathbf{\text{x}}\right),\mathbf{\text{x}}\in R_i^n\hfill \\ \hfill ⋮\hfill \\ \hfill e_l\left(\mathbf{\text{x}}\right),\mathbf{\text{x}}\in R_l^n\hfill \end{array}\right.$ (3)

where e_i(x) is an analytic set of functions within $R_i^n$. It should be noted that at a specific time instance the state is described by a given realization, corresponding to a single branch of equation (3).

Very often the study of the behavior of dynamic systems results in a discretized description of the problem where nodes are connected with each other with components, as can be seen in the following Figure 1. Often the nodes correspond to point masses with the components corresponding to elements that resist the relative displacements and velocities of the masses, or in the case of Finite Element analysis the nodes are the finite element nodes and the components correspond to the finite elements (Zienkiewicz and Taylor, 2005).

FIGURE 1

Figure 1. Systems described as nodal points connected with components (A) Discrete masses connected with components providing resistance to motion (B) nodes connected with finite elements.

The investigated system may be expressed as a combination of individual components C_j, j = 1 … N_c, where N_c is the overall number of the components. For each component C_j a subset, ${\mathbf{\text{x}}}_{C_j}=\left[{\mathbf{\text{x}}}_{C_{j}t},{\mathrm{\theta }}_{C_j}\right]$, of the overall states of the system x, acts as input. The outputs of C_j, ${\mathbf{\text{P}}}_{C_j}$, are connected to the inputs according to a set of equations of the form:

${\dot{\mathbf{\text{P}}}}_{C_j}=E_{C_j}\left({\mathbf{\text{x}}}_{C_j},{\dot{\mathbf{\text{x}}}}_{C_j}\right)$ (4)

where $E_{C_j}$ are generally non-smooth equations, which can further be separated into smooth branches, e.g., for the k^th branch of model C_j:

${\dot{\mathbf{\text{P}}}}_{C_j}=e_{C_j^k}\left({\mathbf{\text{x}}}_{C_j},{\dot{\mathbf{\text{x}}}}_{C_j}\right)$ (5)

where $e_{C_j^k}\left({\mathbf{\text{x}}}_{C_j}\right)$ is an analytic set of functions. As the system evolves dynamically over time, it is expected to shift between the individual branches of a component C_j. This transition between branches will be referred to as a dynamic event, and the corresponding time instance as the time of the event. A set of transition equations $g_{C_j^{k\to l}}\left(\mathbf{\text{x}}\right)=0$ describe the transition from branch $C_j^k$ to $C_j^l$. Having determined the active smooth branch of a component the smooth branch of the overall system from equation (3) can be easily chosen.

2.1. Observability of Non-Smooth Dynamical Systems

The observability of non-smooth systems and the points of differentiation to their smooth counterparts have been discussed in Chatzis et al. (2014, 2017). It should be noted that the notions of observability and identifiability used in these papers and in this work refer to the ability to distinguish the states and parameters from their neighbors at a specific time instance. The method proposed in Chatzis et al. (2014) relies on the study of the observability of each of the smooth subsystems of equation (3), resulting into a characterization of each system branch as either observable, when all associated states are observable, and hence the parameters are identifiable, or as unobservable, when not all states are observable, and hence not all parameters are necessarily identifiable. In general, the separation of an analytic system’s states into an observable and an unobservable set requires a non-linear transformation (Persis and Isidori, 2000). However, for the systems examined herein, it is further assumed that for each of the subsystems i of equation (3), we can further separate the state vector x into its observable and a minimum number of unobservable components, denoted as x^oi and x^ui, in a straightforward manner.

If the union of the observable states from all branches is a strict subset of the state vector $\mathbf{\text{x}}({\cup }_{i=1}^l{\mathbf{\text{x}}}^{\mathit{oi}}\subset \mathbf{\text{x}})$, i.e., does not contain at least one of the components of x, then it may be concluded that these excluded states are unobservable and may not be adequately tracked via a System Identification algorithm. If on the other hand, the union of the observable components results in the state vector x, ${\cup }_{i=1}^l{\mathbf{\text{x}}}^{\mathit{oi}}=\mathbf{\text{x}}$, then each component of the state vector x could potentially be identified within the corresponding smooth branch within which it is observable. Hence, if the response of the system includes at least one branch for which a parameter is identifiable, then a system identification algorithm could potentially succeed in identifying the value of that parameter. In this paper, the latter case of systems is studied, i.e., systems for which the parameters of the model may be inferred via an appropriate system identification method.

For a component C_j whose smooth branch k is defined by equations (5), it is of further interest to proceed in a observability analysis where it is assumed that all of the dynamic states, and their derivatives are measured inputs and P_c are the measured outputs, i.e., $\mathbf{\text{u}}=\left[{\mathbf{\text{x}}}_{C_{j}t},{\dot{\mathbf{\text{x}}}}_{C_{j}t}\right]$ and $y=y_I=\left[{\mathbf{\text{P}}}_{C_j}\right]$. By studying the observability of this system the parameters ${\mathrm{\theta }}_{C_j}$ can be separated into identifiable and a minimum set of unidentifiable parameters $\left[{\mathit{\theta }}_{C_j}^o,{\mathit{\theta }}_{C_j}^u\right]{|}_{y_I}$. Equation (5) can be expressed only in terms of $\left[{{\mathbf{\text{x}}}_{C_{j}t},{\dot{\mathbf{\text{x}}}}_{C_{j}t},\mathrm{\theta }}_{C_j}^o\right]$:

${\dot{\mathbf{\text{P}}}}_{C_j}=e_{C_j^k}\left({\mathbf{\text{x}}}_{C_{j}t},{\dot{\mathbf{\text{x}}}}_{C_{j}t},{\mathit{\theta }}_{C_j}^o\right)$ (6)

Due to the absence of ${\mathit{\theta }}_{C_j}^u{|}_{y_I}$ from equation (6), for any measurement setup which does not directly involve measurement of ${\mathit{\theta }}_{C_j}^u{|}_{y_I}$, those parameters directly contribute to the unidentifiable states x^ui of the corresponding smooth branch of equation (3). If a parameter is shared between different components it will be contributing to the unidentifiable x^ui only if it belongs to ${\mathit{\theta }}_{C_j}^u{|}_{y_I}$ for all of them. It should, however, be noted that whether $\left[{\mathbf{\text{x}}}_{C_{j}t},{\mathit{\theta }}_{C_j}^{\mathbf{\text{o}}}{|}_{{\mathit{\text{y}}}_{\mathit{\text{I}}}}\right]$ contribute to x^ui depends on the observability of the system under the actual measurement setup used.

Hence, this component analysis may often pinpoint part of the unidentifiable parameters, although it ought to further be paired with an observability analysis, as discussed in Chatzis et al. (2014).

3. Unscented Kalman Filter

The UKF simulates non-linear systems by approximating the state as a Gaussian random variable (GRV), represented by a set of carefully chosen deterministic points known as the Sigma Points. This section only provides a basic overview of the filter equations; more details can be found in Julier and Uhlmann (1997) and Wan and Van Der Merwe, 2000 and previous work of the authors Chatzi and Smyth (2009), Chatzi et al. (2010), and Chatzis et al. (2015).

Consider the general dynamical system described by the following equations (7).

${\mathbf{\text{x}}}_k=f\left({\mathbf{\text{x}}}_{k-1},{\mathbf{\text{u}}}_{k-1}\right)+{\mathbf{\text{w}}}_{k-1},{\mathbf{\text{y}}}_k=h\left({\mathbf{\text{x}}}_k,{\mathbf{\text{u}}}_k\right)+{\mathbf{\text{v}}}_k$ (7)

where w_k is the process noise and v_k is the observation noise, both of which are considered to be white Gaussian noise processes of covariance matrices Q_k and P_k, respectively.

Given the state vector at step k − 1 and assuming that this has a mean value of ${\widehat{\mathit{\text{x}}}}_{k-1}$ and covariance P_k−1, the statistics of x_k can be calculated by using the Unscented Transformation, or in other words by computing the set of sigma points ${\mathit{\chi }}_k^i$ with associated weights W_i. The steps of the method are summarized in Table 1.

TABLE 1

Table 1. The steps of the UKF algorithm.

As inferred by the steps outlined above, the UKF algorithm does not discern between observable/unobservable states and identifiable/unidentifiable parameters. The overall convergence of the method is ensured only when a parameter converges faster during identifiable time steps, than it diverges during unidentifiable steps.

4. Discontinuous Unscented Kalman Filter DUKF

As stated in Section 2.1., during a specific time instance where the system lies in a specific smooth branch i, only part x^oi of the state vector may be observable. Therefore, the UKF algorithm is expected to converge only for that observable part x^oi. The predictions furnished during this interval by the UKF for the unobservable part x^ui, which in this work is assumed to be the unidentifiable parameters, are non-optimal and it is also quite likely that during these time intervals the values of x^ui may very well diverge from the real solutions. In fact, their resulting estimates are expected to be inferior to the initial estimates of these parameters in the initiation of the interval. Hence, during such intervals it is argued that the optimal choice would be to update only the observable part of the state.

To introduce the computational part of the D– modification, a row switching transformation matrix T_i is defined such that:

${\mathbf{\text{T}}}_i\mathbf{\text{x}}=\left\{\begin{array}{c}\hfill {\mathbf{\text{x}}}^{\mathit{oi}}\hfill \\ \hfill {\mathbf{\text{x}}}^{\mathit{ui}}\hfill \end{array}\right\}=\mathbf{\text{x}}\prime$ (8)

in other words, pre-multiplying x with T_i results in a rearranged vector x′ where the first n_oi components are observable and the remaining n_ui components are unidentifiable. As the order among the observable components, and likewise for the corresponding order among the unidentifiable parameters, is not of importance any of the T_i matrices that satisfy equation (8) may be chosen. In all cases, those are by definition Boolean matrices containing only one non-zero element per row which is further equal to one, satisfying the property ${\mathbf{\text{T}}}_{\mathbf{\text{i}}}^{-1}={\mathbf{\text{T}}}_{\mathbf{\text{i}}}^{\mathbf{\text{T}}}$. Any vector v and matrix A whose rows, and columns for the latter, correspond to the elements of x may be brought to the order of x′ with the following operations:

$\begin{array}{ll}\hfill {\mathbf{\text{T}}}_i\mathbf{\text{v}}& =\left\{\begin{array}{c}{\mathbf{\text{v}}}^{\mathit{oi}}\\ {\mathbf{\text{v}}}^{\mathit{ui}}\end{array}\right\}{\mathbf{\text{T}}}_i\mathbf{\text{A}}{\mathbf{\text{T}}}_i^T=\left[\sout{\begin{array}{c}\begin{array}{cc}& \\ {\mathbf{\text{A}}}^{\mathit{oo}}& \end{array}\\ \begin{array}{cc}{\mathbf{\text{A}}}^{\mathit{uo}}& \\ & \end{array}\end{array}|\begin{array}{c}\begin{array}{cc}& \\ {\left({\mathbf{\text{A}}}^{\mathit{uo}}\right)}^T& \end{array}\\ \begin{array}{c}{\mathbf{\text{A}}}^{\mathit{uu}}\\ \end{array}\end{array}}\right]\end{array}$ (9)

while for a n × m matrix B, whose rows only correspond to the order of the elements in x, the following operating reorders its elements to the order of x′:

${\mathbf{\text{T}}}_i\mathbf{\text{B}}=\left[\begin{array}{c}\hfill {\mathbf{\text{B}}}^o\hfill \\ \hfill {\mathbf{\text{B}}}^u\hfill \end{array}\right]$ (10)

It is now straightforward to separate a vector or matrix to the observable ^o, unobservable ^u, and cross ^uo components. DUKF follows the steps of the UKF algorithm for steps 1–5, as shown in Table 1. The DUKF structure in nonetheless differentiated from the standard UKF steps as follows: DUKF updates the observable components of the estimates of the mean vector and covariance matrix during the Kalman updating step: ${\widehat{\mathit{\text{x}}}}_{k|k-1}\to {\widehat{\mathit{\text{x}}}}_{k|k}$ and ${\mathbf{\text{P}}}_{k|k-1}\to {\mathbf{\text{P}}}_{k|k}$, using an appropriate Kalman gain matrix defined based on the observable components. The unobservable parts are retained invariant, while the cross terms of the covariance are updated using the Schmidt-Kalman Filter (Schmidt, 1966; Novoselov et al., 2005). The steps of the DUKF algorithm are summarized in detail in Table 2.

TABLE 2

Table 2. The steps of the DUKF algorithm.

The extra steps entailed by the D– modification in comparison to the standard UKF are simple operations involving multiplications with the transformation matrix T_i. Hence, the extra computations required are of minimal cost, while in fact the computational cost of the DUKF results similar to, or lower than, that of the UKF. This may be attributed to the fact that the Kalman updating steps result in multiplications of matrices that are of lower dimension to those of the original method. An additional advantage of the method against the Discontinuous Extended Kalman Filter, DEKF, previously introduced by the authors, lies in that it does not require the detection of events and hence there are no constraints on the algorithm used for the time updating step. This implies that the method can be directly paired with any existing dynamic or finite element software for the time updating step.

4.1. Estimating the Active Smooth Branch

An important step of the DUKF algorithm is related to separating the states to observable and unobservable as indicated in Table 2. This is straightforward to do once the smooth branch the system lies in is known. As discussed earlier, this may more easily be constructed by choosing the active smooth branch for each component C_j of equation (5). If the true value of the states is known that can be done by evaluating the values of the set of functions $g_{C_j}\left(\mathbf{\text{x}}\right)$ which result to the transition equations between branches, $g_{C_j^{k\to l}}\left(\mathbf{\text{x}}\right)=0$. In this paper, this branch has to be estimated by evaluating a related set of functions ${\widehat{g}}_{C_j}\left(\mathbf{\text{x}}\right)$ at $\mathbf{\text{x}}={\widehat{\mathit{\text{x}}}}_{k|k-1}$.

1. The non-smoothness is a result of a non-differentiable function in the state-space equations. In that case, ${\widehat{g}}_{C_j}=g_{C_j}$.

2. The non-smoothness is a result of an inequality constraint equation $g_{C_j}\left(\mathbf{\text{x}}\right)\le \mathbf{\text{0}}$. As a result, at least one of the branches lies entirely within the space defined by the constraint equation $g_{C_j}\left(\mathbf{\text{x}}\right)=\mathbf{\text{0}}$. As for any sigma point i, $g\left({\mathit{\chi }}_{k|k-1}^i\right)<0$ or $g\left({\mathit{\chi }}_{k|k-1}^i\right)=0$, their weighted sum is likely to satisfy: $g_{C_j}\left({\widehat{\mathit{\text{x}}}}_{k|k-1}\right)<0$. In fact the previous could be observed even if $g\left({\mathit{\chi }}_{k|k-1}^i\right)=0,\forall i$. This creates situations where all sigma points may lie in the smooth branch define by the equality constraint, yet the estimated active smooth branch is different. For this reason, ${\widehat{g}}_{C_j}\ne g_{C_j}$ and a different estimator has to be used based on the physics of the studied problem. Such an example is presented in Section 5.2.

5. Applications

5.1. Non-Linear Hysteretic Bouc–Wen Model

In this example, the hysteretic system illustrated in Figure 2 comprising a Bouc–Wen type spring of mass-normalized stiffness k and linear damping c is examined.

FIGURE 2

Figure 2. Mass on a Bouc–Wen Spring.

The relative displacement x of the body with respect to the ground is considered as the measured quantity. The observability of this system was examined in Chatzis et al. (2014). The equations of motion are formulated as:

$\begin{array}{ccc}\hfill & \ddot{x}+kr+c\dot{x}=-{\ddot{x}}_g\hfill & \hfill \\ \hfill & \dot{r}=\dot{x}-\beta \left|\dot{x}\right|{\left|r\right|}^{\mathrm{\nu }-1}r-\mathrm{\gamma }\dot{x}{\left|r\right|}^{\mathrm{\nu }}\hfill \end{array}$ (11)

where k is the stiffness of the spring, c the damping coefficient, and β, γ, and ν are the parameters of the Bouc–Wen model. The term $\dot{r}$ can be rewritten as $\dot{r}=\dot{x}-{\dot{x}}_s$, where x_s is the displacement of the slider and ${\dot{x}}_s=\beta \left|\dot{x}\right|{\left|r\right|}^{\mathrm{\nu }-1}r-\mathrm{\gamma }\dot{x}{\left|r\right|}^{\mathrm{\nu }}$. Hence, r can be thought of as the displacement of the elastic spring. As stated in that paper the dynamic equations of motion of the system can be separated into four smooth branches:

$\begin{array}{ll}& \left(A\right):\dot{r}=\dot{x}-\beta \dot{x}{r}^{\mathrm{\nu }}-\mathrm{\gamma }\dot{x}{r}^{\mathrm{\nu }},\text{for}\dot{x}>0\&r>0\\ & \left(B\right):\dot{r}=\dot{x}+\beta \dot{x}{r}^{\mathrm{\nu }}-\mathrm{\gamma }\dot{x}{r}^{\mathrm{\nu }},\text{for}\dot{x}<0\&r>0\\ & \left(C\right):\dot{r}=\dot{x}+\beta \dot{x}{\left(-r\right)}^{\mathrm{\nu }}-\mathrm{\gamma }\dot{x}{\left(-r\right)}^{\mathrm{\nu }},\text{for}\dot{x}>0\&r<0\\ & \left(D\right):\dot{r}=\dot{x}-\beta \dot{x}{\left(-r\right)}^{\mathrm{\nu }}-\mathrm{\gamma }\dot{x}{\left(-r\right)}^{\mathrm{\nu }},\text{for}\dot{x}<0\&r<0\end{array}$ (12)

within these branches the system is not fully observable, but may be rewritten in the form:

$\begin{array}{ll}& \left(A\right):\dot{r}=\dot{x}-{\mathrm{\Delta }}_1\dot{x}{r}^{\mathrm{\nu }},\text{for}\dot{x}>0\&r>0\\ & \left(B\right):\dot{r}=\dot{x}+{\mathrm{\Delta }}_2\dot{x}{r}^{\mathrm{\nu }},\text{for}\dot{x}<0\&r>0\\ & \left(C\right):\dot{r}=\dot{x}+{\mathrm{\Delta }}_2\dot{x}{\left(-r\right)}^{\mathrm{\nu }},\text{for}\dot{x}>0\&r<0\\ & \left(D\right):\dot{r}=\dot{x}-{\mathrm{\Delta }}_1\dot{x}{\left(-r\right)}^{\mathrm{\nu }},\text{for}\dot{x}<0\&r<0\end{array}$ (13)

where Δ₁ = β + γ and Δ₂ = β − γ. The augmented state vector is hence defined as: $\left[x,\dot{x},r,k,c,\mathrm{\nu },{\mathrm{\Delta }}_1,{\mathrm{\Delta }}_2\right]$. In this new representation, within each branch all of the states $\left(x,\dot{x},r,k,c,\mathrm{\nu }\right)$ are observable while only one of the parameters Δ₁ and Δ₂ is identifiable depending on the $\mathit{sign}\left(\dot{x}r\right)$. When $\dot{x}r\ge 0$, i.e., with branches (A, D) Δ₁ is identifiable, while when $\dot{x}r<0$, i.e., within the branches (B, C) Δ₂ is identifiable.

The previous result can also be demonstrated in terms of the Bouc–Wen spring that can be considered to be the non-smooth component C₁, for which equations (13) correspond to the form of equation (5) with ${\mathbf{\text{P}}}_{C_1}=r$ and ${\mathbf{\text{x}}}_{C_{1}t}=\left[\dot{x}\right]$. For such a component, ${\mathrm{\theta }}_{C_1}^o={\mathrm{\Delta }}_1,{\mathrm{\theta }}_{C_1}^u={\mathrm{\Delta }}_2$, when $\dot{x}r\ge 0$ and ${\mathrm{\theta }}_{C_1}^o={\mathrm{\Delta }}_2,{\mathrm{\theta }}_{C_1}^u={\mathrm{\Delta }}_2$, when $\dot{x}r<0$. For the given measurement setup used, the observability analysis on the overall system shows that there are no further unidentifiable parameters or unobservable states. To complete the description of the method the transformation matrices are defined:

1. $\dot{x}r\ge 0:$

${\mathbf{\text{T}}}_1={\mathbb{I}}_{8\times 8}$ (14)

2. else:

${\text{T}}_2=\left[\sout{\begin{array}{c}\begin{array}{cc}& \\ {\mathbb{I}}_{6\times 6}& \end{array}\\ \begin{array}{cc}{\left\{0\right\}}_{1\times 6}& \\ & \end{array}\end{array}|\begin{array}{c}\begin{array}{cc}& \\ {\left\{0\right\}}_{6\times 1}& \end{array}\\ \begin{array}{cc}0& 1\\ 1& 0\end{array}\end{array}}\right]$ (15)

A system with mass-normalized stiffness and damping terms $k=1000\frac{1}{\mathit{se}{c}^2}$ and $c=2\ast \sqrt{k}\ast 5\%\frac{1}{\mathit{sec}}$, respectively, and Bouc–Wen parameters ν = 2, Δ₁ = 6000, Δ₂ = 2000 initially at rest is subjected to the input ground motion shown in Figure 3. The measured signal is assumed to be the displacement of the system x.

FIGURE 3

Figure 3. (A) Ground acceleration (B) relative displacement.

5.1.1. The Effect of Noise

For the parametric analysis that follows, an initial guess is herein assumed as k₀ = k, c₀ = c, ν = 3, Δ₁/2000 = 2, and Δ₂/2000 = 2. Initially, the effect of different realizations of noise vectors to the convergence of the algorithms will be studied. To that end 1000 different sets of random process and noise vectors are generated, corresponding to a noise-to-signal RMS ratio of 5%. The noisy inputs and outputs are then used in the two methods, DUKF and UKF, which assume corresponding covariance matrices for the process and measurement noise, and the mean error of the final estimates for the Bouc–Wen parameters is calculated. The Cumulative Distribution of the mean BW parameter errors is shown in Figure 4A. Subsequently, the effect of different levels of noise-to-signal RMS ratio for the process noise and the assumed values in the algorithm is investigated. To that end, values of RMS ratios in the range [1%, 7%] with an increment of 1%, where different values are used for the signals that contaminate the input and measurement vector and the assumed covariances of the measurement and process noise used in the models. That creates a total of 7⁴ cases that are examined, for which the mean BW parameter error is calculated, using both the conventional and proposed method, and the results are presented in 4B.

FIGURE 4

Figure 4. Predictions of the two EKF models for the corresponding parameters of the Bouc–Wen model. (A) 5% noise RMS ratio and corresponding assumed covariances. (B) Varying the noise RMS ratio and assumed covariances.

As observed in Figure 4A, DUKF performs superior to the UKF for a given noise-to-signal RMS ratio and is less affected by the exact realizations of the noise vectors. For the DUKF approximately 80% of the cases result into mean parameter error less than 20 versus 60% for the UKF. As can further be deduced from Figure 4B, the same qualitative comparison for the two methods is observed even in the case where diverse noise-to-signal RMS ratios are adopted, while a mismatch is noted between the actual and assumed (in the model) covariances of the process and measurement noise. The use of the DUKF allows for larger discrepancies between the assumed covariances of the noises and their real properties.

5.1.2. The Effect of the Initial Estimates

The effect of different initial estimates for the Bouc–Wen parameters on the convergence is explored in the following Figure 5. To that end, the initial estimates used for Δ₁ and Δ₂ vary in the range between [1, 7] * 2000, while ν is varied in the range: [1.8, 3]. The mean relative error of the BW parameters is calculated for each case and a color is assigned depending on the value of that error. The error color-bar is shown in Figure 5 corresponding to mean errors from 0 to 100%.

FIGURE 5

Figure 5. Mean relative error color-plots for the two methods when different initial values are used for the Bouc–Wen parameters.

As observed in Figure 5, the relative error is lower for the DUKF as compared to the UKF for a wide range of initial estimates of the parameters. Essentially the method is more forgiving in terms of the proximity of the initial estimate to the real value, which offers an important advantage as often the initial estimates are not close to the final value.

5.1.3. Non-Smoothness and Dimensionality

The previous Bouc–Wen spring will be extended to 4 masses connected with Bouc–Wen springs. Each mass is described by a displacement x_i relative to the ground. The non-linear springs are defined by their stiffness k_i and the Bouc–Wen parameters ${\mathrm{\Delta }}_{1_i}$, ${\mathrm{\Delta }}_{2_i}$, and ν_i and linear dampers with coefficients c_i, i = 1, …, 4 as shown in the following Figure 6.

FIGURE 6

Figure 6. 4 masses system with Bouc–Wen springs and linear dampers.

The state-space equations of the system may be assembled after noting that the equation of evolution of the elastic displacement of spring i > 1 becomes:

$\begin{array}{ll}\hfill {\dot{r}}_i& ={\dot{x}}_i-{\dot{x}}_{i-1}-{\mathrm{\Delta }}_{1_i}\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)|r_i{|}^{{\mathrm{\nu }}_i},\text{for}\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)r_i\ge 0\\ \hfill {\dot{r}}_i& ={\dot{x}}_i-{\dot{x}}_{i-1}+{\mathrm{\Delta }}_{2_i}\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)|r_i{|}^{{\mathrm{\nu }}_i},\text{for}\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)r_i<0\end{array}$ (16)

If this system is excited by a ground acceleration ${\ddot{x}}_g$ and the displacements of the four masses [x₁, x₂, x₃, x₄] are measured then it may be demonstrated that for each of the four Bouc–Wen components of the system, C₁, …, C₄, the parameters ${\mathrm{\Delta }}_{1_i}$ and ${\mathrm{\Delta }}_{2_i}$ become unidentifiable when $\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)r_i<0$ or $\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)r_i\ge 0$, respectively. This occurs from the identifiability properties of each Bouc–Wen component, C_i, after noting that ${\mathbf{\text{P}}}_{C_i}=r_i$ and ${\mathbf{\text{x}}}_{C_{i}t}=\left[{\dot{x}}_i-{\dot{x}}_{i-1}\right]$. When $\left({\dot{x}}_i-{\dot{x}}_{i-1}\right)\ge 0$, ${\mathrm{\theta }}_{C_i}^o={\mathrm{\Delta }}_{1_i}$ and ${\mathrm{\theta }}_{C_i}^u={\mathrm{\Delta }}_{2_i}$, else ${\mathrm{\theta }}_{C_i}^o={\mathrm{\Delta }}_{2_i}$ and ${\mathrm{\theta }}_{C_i}^u={\mathrm{\Delta }}_{1_i}$. The remaining parameters and dynamic states are identifiable and observable, respectively, according to the observability analysis of the overall system. The overall transformation matrix T_i can hence be assembled at any time instance.

A system with parameters k_i = [1000, 900, 800, 700] [1/s²], $c_i=2\sqrt{k_i}5∕100$ [1/s], ${\mathrm{\Delta }}_{1_i}=\left[6,7,8,9\right]$ 2000, ${\mathrm{\Delta }}_{1_i}=\left[6,7,8,9\right]$ 2000, ${\mathrm{\Delta }}_{1_i}=\left[1,1,2,6\right]$ 1000, and ν_i = [2, 2, 2, 2] is subjected to the time history of Figure 3A. The obtained displacements are shown in the following Figure 7.

FIGURE 7

Figure 7. Displacements of the four masses of the system.

The measurements are contaminated with noise of noise-to-signal RMS ratio of 1%. The corresponding covariances are assumed in both models for the process and measurement noises. The initial estimates used for both models are k_i = 1000, $c_i=2\sqrt{1000}10∕100$, ν_i = 2.5, ${\mathrm{\Delta }}_{1_i}={\mathrm{\Delta }}_{2_i}=3$ for i = 1, ⋅, 4. The results of the identification using the UKF and DUKF are shown in the following Figures 8 and 9.

FIGURE 8

Figure 8. Spring and damping parameters identified by the UKF and DUKF.

FIGURE 9

Figure 9. Bouc–Wen parameters identified by the UKF and DUKF.

As observed in Figure 8, both methods provide fairly good estimates of the elastic parameters of the system k_i, c_i. However, when it comes to the non-linear (Bouc–Wen) parameters, the DUKF provides a substantially improved estimated versus the UKF. This can be seen by the fact that while for DUKF the final ratio of the estimated over real values for the parameters is close to unity, for the UKF this ratio substantially deviates from unity indicating large estimation errors. This can be attributed to the multiple, 4 in this case, unidentifiable parameters at any time window. As those unidentifiable parameters are increased, it is more likely that the estimates of the system overall diverge.

It should hence be noted that non-smooth high dimensional systems suffer from the effects of dimensionality, but also from the additional effect of the increased number of unidentifiable parameters and hence sources of divergence. The former can be improved using techniques applied to smooth systems for dimensionality (Olivier and Smyth, 2017b), while the latter is treated through the D– modification suggested here. It should be noted that the two treatments, which aim at tackling different problems, can be combined.

5.2. 2DOF Elasto-Plastic System

In this example, the behavior of a shear system of two masses with displacements x₁ and x₂ connected to each other and the ground by means of linear damping elements of normalized damping over mass c₁, c₂ and elastoplastic springs of normalized over mass stiffness k₁, k₂ and yield force $F_{y_1}$ and $F_{y_2}$ as shown in Figure 10 is studied.

FIGURE 10

Figure 10. System studied. (A) System of two masses connected with elastoplastic springs and linear dampers. (B) Behavior of elastoplastic spring.

The equations of motion describing the system when subjected to a ground acceleration ${\ddot{x}}_g$ become:

$\begin{array}{ll}\hfill {\ddot{x}}_1+\left(c_1+c_2\right){\dot{x}}_1-c_2{\dot{x}}_2+k_1{x}_{{\mathit{el}}_1}-k_2{x}_{{\mathit{el}}_2}& =-{\ddot{x}}_g\\ \hfill {\ddot{x}}_2+\left(c_2\right){\dot{x}}_2-c_2{\dot{x}}_1+k_2{x}_{{\mathit{el}}_2}& =-{\ddot{x}}_g\end{array}$ (17)

where $x_{e{l}_i}$ is the elastic elongation of the elastoplastic spring i whose evolution over time is defined as:

${\dot{x}}_{{\mathit{el}}_1}={\dot{x}}_1,\text{in the elastic branch}{\dot{x}}_{{\mathit{el}}_1}=0,\text{in the plastic branch}$ (18)

${\dot{x}}_{{\mathit{el}}_2}={\dot{x}}_2-{\dot{x}}_1,\text{in the elastic branch}{\dot{x}}_{{\mathit{el}}_2}=0,\text{in the plastic branch}$ (19)

The following equations define the transition conditions between the elastic and plastic branches for spring i:

$\begin{array}{ll}\hfill \parallel k_i{x}_{{\mathit{el}}_i}\parallel & =F_{y_i},\text{elastic}\to \text{plastic}\\ \hfill {\dot{x}}_{{\mathit{el}}_i}& =0,\text{plastic}\to \text{elastic}\end{array}$ (20)

The previous transition equations ensure that the force of elastoplastic spring i always satisfies the condition: $\parallel k_i{x}_{e{l}_i}\parallel \le F_{y_i}$.

For the purpose of identification the augmented state vector will include $x_1,x_2,{\dot{x}}_1,{\dot{x}}_2,\mathit{kxe}{l}_1,\mathit{kxe}{l}_2,c_1,c_2,k_1,k_2,F_{y_1},F_{y_2}$. The dynamic states $\mathit{kxe}{l}_i$ correspond to the product $k_i{x}_{e{l}_i}$ and are used instead of the elastic displacements as it allows separating the states into observable and unidentifiable within all branches without having to use a non-linear transformation (Chatzis et al., 2017). In terms of the identifiability of the system it can easily be shown as in Chatzis et al. (2017), that all the dynamic states together with c₁ and c₂ are always identifiable. However, only one of $\left[k_i,F_{y_i}\right]$ is identifiable depending on whether spring i lies in an elastic or plastic branch at that specific time instance. This follows after studying the identifiability of any of the two elastoplastic spring components, C_i where ${\mathbf{\text{P}}}_{C_i}=\mathit{kxe}{l}_i$ and $x_{C_{1}t}=\left[{\dot{x}}_i-{\dot{x}}_{i-1}\right]$. Then the equation of component C_i becomes:

$\begin{array}{ll}\hfill {\dot{\mathbf{\text{P}}}}_{C_i}& =k_i\left({\dot{x}}_i-{\dot{x}}_{i-1}\right),\text{in the elastic branch}\\ \hfill {\dot{\mathbf{\text{P}}}}_{C_i}& =0,\text{in the plastic branch}\end{array}$ (21)

and as a result ${\mathrm{\theta }}_{C_i}^u=F_{y_i}$ in the elastic branch and ${\mathrm{\theta }}_{C_i}^u=k_i$ in the plastic branch. $F_{y_i}$ becomes identifiable when component C_i enters the plastic branch. This is because the plasticity constraint, when activated, effectively results into an additional measurement equation: $F_{y_i}=\parallel \mathit{kxe}{l}_i\parallel$.

Hence, the identifiability of the system requires estimation of whether spring i is in an elastic or plastic branch. While this discussion is obvious for the real system through use of equation (20), it requires some careful consideration when applied to the systems estimated by the DUKF. As each of the sigma points are bound by the inequality constraint of equation (20), it is likely that their mean would satisfy the condition $\parallel {k\widehat{x}\mathit{el}}_i\parallel <F_{y_i}$ even if the majority of the sigma points are satisfying the equality, and are hence in the plastic branch. The problem occurs due to the fact that the condition in equation (20) indicates of whether the spring is elastic or plastic, but cannot quantify “how” elastic or plastic the response of the system is. A means of indicating the tendency of the system to behave in an elastic or plastic manner, suggested in this paper, may be attained via comparison of the estimated mean velocities of the elastic and plastic elongation of the springs.

To such an end, the inequality of equation (20) is used to deem of whether each sigma point lies in an elastic or plastic branch using the estimated values for the springs prior to applying the measurement update (i.e., ${\mathrm{\chi }}_{k|k-1}^i$ for sigma point i). Then, the elastic and plastic velocity of spring i for sigma point j, ${\dot{x}}_{e{l}_j}^j$, ${\dot{x}}_{p{l}_j}^j$, are:

$\begin{array}{ll}\hfill {\dot{x}}_{{\mathit{el}}_i}^j& ={\dot{x}}_{s_i}^j,\text{if spring}i\text{is elastic}\\ \hfill {\dot{x}}_{{\mathit{el}}_i}^j& =0,\text{if spring}i\text{is plastic}\end{array}$ (22)

where ${\dot{x}}_{s_i}^j$ is the total velocity of spring i for sigma point j, ${\dot{x}}_{s_1}^j={\dot{x}}_1^j$ and ${\dot{x}}_{s_2}^j={\dot{x}}_2^j-{\dot{x}}_1^j$, then ${\dot{x}}_{p{l}_i}^j={\dot{x}}_{s_i}^j-{\dot{x}}_{e{l}_i}^j$. The mean estimates of the two velocities ${\widehat{\dot{x}}}_{e{l}_i}$ and ${\widehat{\dot{x}}}_{p{l}_i}$ can be calculated using the fourth step of the DUKF algorithm. Hence, the following criterion is used to deem of whether the system is behaving elastically or plastically:

$\begin{array}{ll}\hfill \parallel {\widehat{\dot{x}}}_{{\mathit{el}}_i}\parallel \ge \parallel {\widehat{\dot{x}}}_{{\mathit{pl}}_i}\parallel & \to \text{the system behaves elastically}\\ \hfill \parallel {\widehat{\dot{x}}}_{{\mathit{el}}_i}\parallel <\parallel {\widehat{\dot{x}}}_{{\mathit{pl}}_i}\parallel & \to \text{the system behaves plastically}\end{array}$ (23)

It is now straightforward to estimate the branch each spring would be in and obtain the corresponding contribution to the transformation matrix T. It should finally be noted that in both the UKF and DUKF algorithms the following constraint is applied to sigma point j if $\parallel \mathit{kxe}{l}_i\parallel >F_{y_i}$ for spring i:

${\mathit{kxel}}_i=F_{y_i}\ast \mathit{sign}\left({\mathit{kxel}}_i\right)$ (24)

where equation (24) is a return mapping scheme. This is what one would follow in the forward dynamics problem, as the value of $F_{y_i}$ would be known. However, in this problem it would also be possible to instead modify the value of $F_{y_i}$ when the plasticity constraint is violated.

A system with properties k₁ = 1000 [1/s²], k₂ = 800 [1/s²], $c_1=2\sqrt{k_1}0.05$ [1/s], $c_2=2\sqrt{k_2}0.05$ [1/s], $F_{y_1}=50$ [m/s²] $F_{y_1}=30$ [m/s²] is subjected to the excitation of Figure 11A. The occurring displacements of both masses are measured as shown in Figure 11B.

FIGURE 11

Figure 11. (A) Ground acceleration and (B) displacements of the masses.

The two springs exhibit an elastoplastic response as shown in the force displacement responses plotted in Figure 12. The maximum total displacements, x₁ and x₂ − x₁ correspond to 1.7 and 2 times the yield displacements, respectively.

FIGURE 12

Figure 12. Force displacement curves for (A) Spring 1 and (B) Spring 2.

The input and measurements are contaminated with white noise signals corresponding to 5% noise-to-signal RMS ratios. As there is a substantial drift of the measured signals their RMS is calculated after these are passed through a high pass filter with a cutoff frequency at 0.5 Hz. The initial estimates for the stiffness and the damping of both springs are given a significant offset, as these are assumed as twice their actual value. The initial estimates of the yield forces, $F_{y_1}$ and $F_{y_2}$, are varied in the following ranges: $F_{y_1}\in \left[5,80\right]$ and $F_{y_2}\in \left[5,60\right]$. This is later used in both the UKF and DUKF, where the assumed covariances of the process and measurement noises match the 5% noise-to-signal RMS ratio. After the algorithms are implemented, the mean relative error of the estimated parameters with respect to the real values is calculated and is plotted in the following Figure 13, where the upper row of figures corresponds to the results of DUKF and the lower to UKF, each column of sub-figures corresponds to the error of the parameter indicated by the title and for each sub-figure the horizontal and vertical axes correspond to different initial estimates of $F_{y_1}$ and $F_{y_2}$. The mean relative error is indicated by the color-bar of the Figure.

FIGURE 13

Figure 13. Mean relative error color-plots for the two methods when varying the initial estimates of ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{1}}}$ and ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{2}}}$.

As observed in Figure 13, DUKF in general results in reduced errors over UKF for a wide range of initial estimates. It should be noted that the method results in a large improvement over the estimates of the stiffness of the two springs k₁ and k₂. This is expected as these parameters become unidentifiable when the corresponding spring is in the plastic branch. Equally there appears to be a clear improvement for the estimates of the plastic forces $F_{y_1}$ and $F_{y_2}$, when the initial estimates are in the range ${F_y}_1\in \left(0,70\right)$ and $F_{y_2}\in \left(5,50\right)$. For values of $F_{y_1}>70$ and $F_{y_2}>50$, DUKF does not change the initial estimate of the corresponding parameter, as the algorithms estimates that the system is always elastic. This can be understood by looking at the following Figure 14 which is plotting the forces seen by elastic springs of stiffness k₁ and k₂ for the real displacements of the system.

FIGURE 14

Figure 14. Force displacement curves for (A) Spring 1 and (B) Spring 2.

In both cases, there are only few points in Figures 14A,B, where, respectively, even the force of a linear spring for the displacements of the system would exceed the values of $F_{y_1}>70$ and $F_{y_2}>50$. As a result, it is reasonable for the DUKF, given the properties of the estimator selected and the way the plasticity constraint is applied, to reach the conclusion that the corresponding spring always remain elastic when the initial estimates of $F_{y_i}$ are in the aforementioned range. While this is disadvantageous in terms of identification of the real value of the parameters, when high initial estimates of $F_{y_1}$ and $F_{y_2}$ are used, the advantage of the DUKF lies in that the algorithm has not changed the estimates of the corresponding covariance terms indicating that these parameters were not identified. Additionally, it appears that even for such cases the DUKF is capable of providing good estimates for the remaining parameters.

In contrast, the UKF would proceed to evolve the initial estimate of $F_{y_1}$ and $F_{y_2}$ in all scenarios. The algorithm may hence reduce the initial estimate of $F_{y_1}$ or $F_{y_2}$ even during periods of unidentifiability and this non-optimal change could happen to result into more favorable estimates for the value of $F_{y_1}$ and $F_{y_2}$. However, it is equally or more probable that the algorithm will change the overall estimates of the parameters to less favorable values during unidentifiable windows, thus resulting in divergence of all parameters. This is depicted in the behavior of the UKF in Figure 13 for initial estimates of $F_{y_1}>70$ and $F_{y_2}>50$ where even when the algorithm happens to do better for the estimates of the corresponding $F_{y_i}$ the remaining parameters behave less optimally. Additionally, in that region the final estimate of the covariance terms corresponding to $F_{y_i}$ are substantially lower than the DUKF even when the algorithm does not converge. The behavior of both algorithms for the case of initial estimates $F_{y_1}=70$, $F_{y_2}=60$ is shown in the following Figures 15 and 16 for the estimated/real values of the parameters versus time and the estimated versus real plastic displacements $x_{p{l}_i}$, respectively.

FIGURE 15

Figure 15. Estimated over real values for the parameters, for the disadvantageous for the DUKF scenario of initial estimates ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{1}}}=\text{70}$, ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{2}}}=\text{60}$.

FIGURE 16

Figure 16. Estimated versus real plastic displacements of the springs for (A) UKF and (B) DUKF, for the disadvantageous for the DUKF scenario of initial estimates ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{1}}}=\text{70}$, ${\mathit{\text{F}}}_{{\mathit{\text{y}}}_{\text{2}}}=\text{60}$.

It should be noted that in Figure 15, the DUKF does not practically alter the initial guess of $F_{y_i}$ as described above. In this simulation the UKF happens to evolve the estimate of $F_{y_1}$ in a coincidentally favorable way, while doing the opposite happens for $F_{y_2}$. However, as observed in Figure 16A, the UKF has difficulty in tracking the plastic displacements for this case. To the contrary, the DUKF is shown in Figure 16B to result into excellent predictions despite the inability to track the values of $F_{y_1}$ and $F_{y_2}$. Finally it should be noted that UKF appears relatively certain for the values of $F_{y_1}$ and $F_{y_2}$, despite the fact that the latter is incorrect, with corresponding covariance terms smaller than 2 × 10⁻³. In contrast, DUKF has not practically changed the covariance from the initial guess which is of the order of 10² indicating that the method is uncertain of its estimation of these two values.

Hence, this investigation leads to a result similar to those of the DEKF in Chatzis et al. (2017) for Elasto-plastic springs: it appears more favorable to use initial estimates of $F_{y_i}$ smaller than the real value of those parameters, or at least smaller than the maximum force seen by an elastic spring for the displacements of the system. Of course while this is not possible a priori, as those values are not known one can be informed by the DUKF of the inability of the algorithm to identify the corresponding value of $F_{y_i}$, due to the lack of change of the parameters and the corresponding large covariance terms. Additionally, the elastic displacements $x_{e{l}_i}$ and total displacements x_i of the system are calculated by the DUKF with high precision allowing to alert the user on the presence of permanent displacements.

6. Discussion and Conclusion

This paper suggests the use of a Discontinuous D– modification for non-smooth systems for modifying the UKF. Non-smooth systems include certain parameters whose identifiability property changes over time. To alleviate the divergence exhibited by standard filtering algorithms during time periods of unidentifiability, the Discontinuous modification suggests retaining such parameters invariant during those intervals. This paper, therefore, introduces a Discontinuous Unscented Kalman Filter DUKF.

The method is implemented as a minimally invasive modification allowing, as the original UKF algorithm, straightforward implementation with any existing software that employ filtering algorithms to update the states of a system over time. The D– modification makes use of a transformation matrix at any time instance based on the estimated active smooth branch of the system and the occurring identifiable and unidentifiable parameters. The proposed modification does not increase the computational cost of the original method; in fact it leads to inversions and operations on matrices of lower dimension.

The examples provided demonstrate the use of the algorithm with two different types of non-smooth systems: systems where the identifiability condition varies between different subspaces of the state vector and systems for which the non-smoothness is a result of an inequality constraint. Several non-smooth systems can be described as combinations of these two cases. The examples illustrate the robustness of the D– modification and the overall improved performance of DUKF versus UKF for problems of increased complexity. Different sources of complexity were studied in terms of their effect such as the noise in the input and measured data, the assumed noise covariances in the models, different initial conditions and dimensionality. The DUKF was shown to improve the estimates provided for all the previous cases, and resulted in a more consistent behavior than that of the standard UKF for non-smooth systems.

This paper together with former work of the authors on the EKF (Chatzis et al., 2017) demonstrate that non-smoothness bears an effect on the convergence of non-linear Kalman Filters and in general for online Bayesian methods, further illustrating that the D– modification is a viable treatment across algorithmic implementations of this class. The D– modification tackles the problems associated to non-smoothness and may be paired with existing modifications aiming at improving the performance of the algorithms for smooth problems, substantially expanding the ability of the occurring algorithms to handle problems of increased complexity.

Author Contributions

Both authors verify to have contributed to this original research paper, which has not been submitted for publication elsewhere.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

MC would like to acknowledge the Marie Curie FP7 Career Integration Grant No. 618359 within the seventh European Union Framework Programme, for the support of this research. EC would like to gratefully acknowledge the support of the Albert Lück Stiftung. The authors would also like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work: http://dx.doi.org/10.5281/zenodo.2255.

References

Alavi, S., Mahdi, A., Jacob, P. E., Payne, S. J., and Howey, D. A. (2015). Structural Identifiability Analysis of Fractional Order Models with Applications in Battery Systems. arXiv preprint arXiv:1511.01402.

Google Scholar

Al-Hussein, A., and Haldar, A. (2015). Novel unscented Kalman filter for health assessment of structural systems with unknown input. J. Eng. Mech. 141, 4015012.

Google Scholar

Astroza, R., Ebrahimian, H., and Conte, J. P. (2017). Batch and Recursive Bayesian Estimation Methods for Nonlinear Structural System Identification. Springer, 341–364.

Google Scholar

Au, S.-K., and Zhang, F.-L. (2016). Fundamental two-stage formulation for Bayesian system identification, part I: general theory. Mech. Syst. Signal Process. 66, 31–42. doi: 10.1016/j.ymssp.2015.04.025

CrossRef Full Text | Google Scholar

Chatzi, E., and Smyth, A. (2009). The unscented Kalman filter and particle filter methods for nonlinear structural system identification with non-collocated heterogeneous sensing. Struct. Control Health Monit. 16, 99–123. doi:10.1002/stc.290

CrossRef Full Text | Google Scholar

Chatzi, E. N., and Smyth, A. W. (2012). Particle filter scheme with mutation for the estimation of time-invariant parameters in structural health monitoring applications. Struct. Control Health Monit. 20, 1081–1095. doi:10.1002/stc.1520

CrossRef Full Text | Google Scholar

Chatzi, E. N., and Smyth, A. W. (2014). Nonlinear System Identification: Particle-Based Methods. Berlin, Heidelberg: Springer, 1–18.

Google Scholar

Chatzi, E. N., Smyth, A. W., and Masri, S. F. (2010). Experimental application of on-line parametric identification for nonlinear hysteretic systems with model uncertainty. J. Struct. Saf. 32, 326–337. doi:10.1016/j.strusafe.2010.03.008

CrossRef Full Text | Google Scholar

Chatzis, M., Chatzi, E., and Smyth, A. (2015). An experimental validation of time domain system identification methods with fusion of heterogeneous data. Earthquake Eng. Struct. Dyn. 44, 523–547. doi:10.1002/eqe.2528

CrossRef Full Text | Google Scholar

Chatzis, M., Chatzi, E., and Triantafyllou, S. (2017). A discontinuous extended Kalman filter for non-smooth dynamic problems. Mech. Syst. Signal Process. 92, 13–29. doi:10.1016/j.ymssp.2017.01.021

CrossRef Full Text | Google Scholar

Chatzis, M. N., Chatzi, E. N., and Smyth, A. W. (2014). On the observability and identifiability of nonlinear structural and mechanical systems. Struct. Control Health Monit. 22, 574–593.

Google Scholar

Chatzis, M. N., and Smyth, A. W. (2012a). Modeling of the 3d rocking problem. Int. J. Nonlinear Mech. 47, 85–98. doi:10.1016/j.ijnonlinmec.2012.02.004

CrossRef Full Text | Google Scholar

Chatzis, M. N., and Smyth, A. W. (2012b). Robust modeling of the rocking problem. J. Eng. Mech. 3, 247–262. doi:10.1061/(ASCE)EM.1943-7889.0000329

CrossRef Full Text | Google Scholar

Chatzis, M. N., and Smyth, A. W. (2013). Three-dimensional dynamics of a rigid body with wheels on a moving base. J. Eng. Mech. 139. doi:10.1061/(ASCE)EM.1943-7889.0000456

CrossRef Full Text | Google Scholar

Ding, Y., Zhao, B. Y., and Wu, B. (2014). Structural system identification with extended Kalman filter and orthogonal decomposition of excitation. Math. Probl. Eng. 2014, 10. doi:10.1155/2014/987694

CrossRef Full Text | Google Scholar

Diop, S., and Fliess, M. (1991). “Nonlinear observability, identifiability, and persistent trajectories,” in Decision and Control, 1991, Proceedings of the 30th IEEE Conference on (Brighton: IEEE), 714–719.

Google Scholar

Ebrahimian, H., Astroza, R., and Conte, J. P. (2015). Extended Kalman filter for material parameter estimation in nonlinear structural finite element models using direct differentiation method. Earthquake Eng. Struct. Dyn. 44, 1495–1522. doi:10.1002/eqe.2532

CrossRef Full Text | Google Scholar

Ebrahimian, H., Astroza, R., Conte, J. P., and de Callafon, R. A. (2017). Nonlinear finite element model updating for damage identification of civil structures using batch Bayesian estimation. Mech. Syst. Signal Process. 84(Part B), 194–222. doi:10.1016/j.ymssp.2016.02.002

CrossRef Full Text | Google Scholar

Eftekhar Azam, S., Ghisi, A., and Mariani, S. (2012). Parallelized sigma-point Kalman filtering for structural dynamics. Comput. Struct. 92-93, 193–205. doi:10.1016/j.compstruc.2011.11.004

CrossRef Full Text | Google Scholar

Farrar, C. R., and Worden, K. (2012). Structural Health Monitoring: A Machine Learning Perspective. Chichester: John Wiley & Sons, Ltd.

Google Scholar

Giannakopoulos, A. E. (1989). The return mapping method for the integration of friction constitutive relations. Comput. Struct. 32, 157–167. doi:10.1016/0045-7949(89)90081-3

CrossRef Full Text | Google Scholar

Greenbaum, R. J., Smyth, A. W., and Chatzis, M. N. (2015). Monocular computer vision method for the experimental study of three-dimensional rocking motion. J. Eng. Mech. 142, 04015062. doi:10.1061/(ASCE)EM.1943-7889.0000972

CrossRef Full Text | Google Scholar

Hermann, R., and Krener, A. (1977). Nonlinear controllability and observability. IEEE Trans. Autom. Control 22, 728–740. doi:10.1109/TAC.1977.1101601

CrossRef Full Text | Google Scholar

Huang, J., McBratney, A. B., Minasny, B., and Triantafilis, J. (2017). Monitoring and modelling soil water dynamics using electromagnetic conductivity imaging and the ensemble Kalman filter. Geoderma 285, 76–93. doi:10.1016/j.geoderma.2016.09.027

CrossRef Full Text | Google Scholar

Julier, S. J., and Uhlmann, J. K. (1997). “A new extension of the Kalman filter to nonlinear systems,” in Proceedings of AeroSense: The 11th Int. Symposium on Aerospace/Defense Sensing, Simulation and Controls (Orlando).

Google Scholar

Kakouris, E., and Triantafyllou, S. (2017). Material point method for crack propagation in anisotropic media: a phase field approach. Arch. Appl. Mech. 1–30. doi:10.1007/s00419-017-1272-7

CrossRef Full Text | Google Scholar

Kalman, R. (1963). Mathematical description of linear dynamical systems. J. Soc. Ind. Appl. Math. A Control 1, 152–192. doi:10.1137/0301010

CrossRef Full Text | Google Scholar

Kumar, A., and Crassidis, J. (2007). “Colored-noise Kalman filter for vibration mitigation of position/attitude estimation systems,” in Guidance, Navigation, and Control and Co-Located Conferences (Hilton Head, South Carolina: American Institute of Aeronautics and Astronautics).

Google Scholar

Liu, P.-T., Li, F., and Xiao, H. (1996). A state decoupling approach to estimate unobservable tracking systems. IEEE J. Oceanic Eng. 21, 256–259. doi:10.1109/48.508156

CrossRef Full Text | Google Scholar

Ljung, L., and Glad, T. (1994). On global identifiability for arbitrary model parametrizations. Automatica 30, 265–276. doi:10.1016/0005-1098(94)90029-9

CrossRef Full Text | Google Scholar

Mariani, S., and Corigliano, A. (2005). Impact induced composite delamination: state and parameter identification via joint and dual extended Kalman filters. Comput. Methods Appl. Mech. Eng. 194, 5242–5272. doi:10.1016/j.cma.2005.01.007

CrossRef Full Text | Google Scholar

Novoselov, R. Y., Herman, S. M., Gadaleta, S. M., and Poore, A. B. (2005). “Mitigating the effects of residual biases with Schmidt-Kalman filtering,” in Information Fusion, 2005 8th International Conference on, Vol. 1 (Sweden: IEEE), 8.

Google Scholar

Olivier, A., and Smyth, A. W. (2017a). On the performance of online parameter estimation algorithms in systems with various identifiability properties. Front. Built Environ. 3:14. doi:10.3389/fbuil.2017.00014

CrossRef Full Text | Google Scholar

Olivier, A., and Smyth, A. W. (2017b). Particle filtering and marginalization for parameter identification in structural systems. Struct. Control Health Monit. 24, e1874. doi:10.1002/stc.1874

CrossRef Full Text | Google Scholar

Omrani, R., Hudson, R. E., and Taciroglu, E. (2013). Parametric identification of nondegrading hysteresis in a laterally and torsionally coupled building using an unscented Kalman filter. J. Eng. Mech. 139. doi:10.1061/(ASCE)EM.1943-7889.0000498

CrossRef Full Text | Google Scholar

Papadimitriou, C., and Papadioti, D.-C. (2013). Component mode synthesis techniques for finite element model updating. Comput. Struct. 126, 15–28. doi:10.1016/j.compstruc.2012.10.018

CrossRef Full Text | Google Scholar

Persis, C. D., and Isidori, A. (2000). On the observability codistributions of a nonlinear system. Syst. Control Lett. 40, 297–304. doi:10.1016/S0167-6911(00)00014-1

CrossRef Full Text | Google Scholar

Schmidt, S. F. (1966). Applications of state space methods to navigation problems. Adv. Control Syst. 3, 293–340. doi:10.1016/B978-1-4831-6716-9.50011-4

CrossRef Full Text | Google Scholar

Smyth, A. W., Masri, S. F., Chassiakos, A. G., and Caughey, T. K. (1999). On-line parametric identification of MDOF nonlinear hysteretic systems. J. Eng. Mech. 125, 133–142. doi:10.1061/(ASCE)0733-9399(1999)125:2(133)

CrossRef Full Text | Google Scholar

Villaverde, A. F., and Banga, J. R. (2017). Structural properties of dynamic systems biology models: identifiability, reachability, and initial conditions. Processes 5, 29. doi:10.3390/pr5020029

CrossRef Full Text | Google Scholar

Villaverde, A. F., Barreiro, A., and Papachristodoulou, A. (2016). Structural identifiability of dynamic systems biology models. PLoS Comput. Biol. 12:e1005153. doi:10.1371/journal.pcbi.1005153

PubMed Abstract | CrossRef Full Text | Google Scholar

Walter, E. (1982). Identifiability of State Space Models. Berlin, Heidelberg: Springer.

Google Scholar

Wan, E., and Van Der Merwe, R. (2000). “The unscented Kalman filter for nonlinear estimation,” in Adaptive Systems for Signal Processing, Communications, and Control Symposium, AS-SPCC (Lake Louise, Alberta, Canada: IEEE), 153–158.

Google Scholar

Worden, K., Farrar, C. R., Haywood, J., and Todd, M. (2008). A review of nonlinear dynamics applications to structural health monitoring. Struct. Control Health Monit. 15, 540–567. doi:10.1002/stc.215

CrossRef Full Text | Google Scholar

Wriggers, P. (1991). Computational Contact Mechanics. Berlin, Heidelberg: Springer.

Google Scholar

Zhang, F.-L., and Au, S.-K. (2016). Fundamental two-stage formulation for Bayesian system identification, part II: application to ambient vibration data. Mech. Syst. Signal Process. 66, 43–61. doi:10.1016/j.ymssp.2016.03.024

CrossRef Full Text | Google Scholar

Zienkiewicz, O. C., and Taylor, R. L. (2005). The Finite Element Method for Solid and Structural Mechanics, 6th Edn. Waltham, MA: Elsevier.

Google Scholar