A Bayesian optimized framework for successful application of unscented Kalman filter in parameter identification of MDOF structures

Abstract

The success of the unscented Kalman filter can be jeopardized if the required initial parameters are not identified carefully. These parameters include the initial guesses and the levels of uncertainty in the target parameters and the process and measurement noise parameters. While a set of appropriate initial target parameters give the unscented Kalman filter a head start, the uncertainty levels and noise parameters set the rate of convergence in the process. Therefore, due to the coupling effect of these parameters, an inclusive approach is desired to maintain the chance of convergence for expensive experimental tests. In this paper, a framework is proposed that, via a virtual emulation prior to the experiment, determines a set of initial conditions to ensure a successful application of the online parameter identification. A Bayesian optimization method is proposed, which considers the level of confidence in the initial guesses for the target parameters to suggest the appropriate noise covariance matrices. The methodology is validated on a five-story shear frame tested on a shake table. The results indicate that, indeed, a trade-off can be made between the robustness of the online updating and the final parameter accuracy.

Introduction

Model updating techniques are widely used in structural health monitoring applications such as damage detection, ¹ force identification, ² response prediction, ³ and structural control of the existing structures. ⁴ Dynamic model updating calibrates numerical models to improve the simulation accuracy of the structural responses under certain loads. The calibration is achieved through the system identification of the structure under study. Although linear system identification methods have been studied and successfully applied to various structures, ⁵ due to the nonlinearities in real-world structures, the range of applicability is limited for these methods. Kalman filter is a member of the Bayes filter family, which is frequently used for state estimation of linear models and Gaussian distributed parameters. However, for nonlinear systems, to avoid obtaining non-Gaussian distributions for the transformed parameters of the state vector, either a linearization technique or a non-parametric method should be adopted.

Various nonlinear system identification techniques have been suggested in the literature, such as the least square estimation, ⁶ ${\mathrm{H}}_{\mathrm{\infty }}$ Filter, ⁷ the extended Kalman filter (EKF), ⁸ the unscented Kalman filter (UKF),^9–12 and sequential Monte Carlo^13,14 methods. Among them, UKF and EKF have become popular due to their practicality and computational efficiencies for Gaussian-distributed parameters.^15–18 EKF was first introduced to extend the Kalman filter applications to nonlinear models by locally linearizing the functions at each step of the process. ¹⁹ EKF has been successfully applied to models with low-to-moderate nonlinearities.^20–22 However, since the linearization is performed by a Taylor expansion, and Jacobian matrices are calculated at each step, EKF may become inefficient for high-dimensional variable estimation and large uncertainties.

The UKF, on the other hand, does not require the linearization of the underlying functions and, instead, propagates the probability distribution functions (PDF) of the parameters in an efficient way. UKF directly transforms a set of weighted sigma points through the nonlinear function and approximates the Gaussian distribution of the state and measurement vectors. Comparative studies have concluded the better robustness of UKF compared to EKF.^10,23 UKF is more robust when noisy measurements are present and is accurate up to the second-order Taylor series expansion, which, in case of Gaussian variables, handles up to the third order for any nonlinearity.^24,25

UKF has been implemented in various system identification problems in civil engineering. The computational efficiency of UKF has made it favorable for online updating methods such as real-time hybrid simulation (RTHS).^26–29 Various numerical studies have been conducted on different aspects of the UKF algorithm. Shao et al. ³⁰ implemented UKF for RTHS of a 3 degree of freedom (DOF) structure. The importance of initial covariance matrices of state and noise vectors, as well as the initial guesses for the structural parameters that are the focus of the system identification (called the target parameters hereafter), was stated. Parametric study of the mentioned variables was conducted for the specific case, and the initial target parameters for the RTHS experiment were obtained from a prior quasi-static experimental testing.

Variants of the UKF method have been suggested to overcome different possible scenarios. Constrained unscented Kalman filter is proposed for enforcing boundaries on the state variables for higher robustness of the algorithm. ³¹ Adaptive UKF is proposed for online identification of sudden changes in stiffness of structures under seismic loads. ³² Xie and Feng ¹⁵ studied the applications of the iterated unscented Kalman filter (IUKF) to highly nonlinear systems subjected to different signal noise levels and concluded that the IUKF is more accurate in highly nonlinear systems and more robust to measurement noise levels. Al-Hussein and Haldar ³³ suggested a substructure concept which estimates the stiffness and Rayleigh damping coefficients for a small part of the structure and generates the excitation using an iterative least-squares method. Responses of the key nodes of the structure are required for this implementation, and the stiffness and damping coefficients are used as the initial state vector for the UKF procedure. The authors later suggested a weighted global iteration (WGI) UKF, which, instead of implementing the UKF to the complete set of responses once, applies the UKF to a fraction of the responses several times and reduces the uncertainty of the state parameters. ³⁴ Although the excitation data is no longer required in this approach, the system responses such as velocity and displacement cannot be obtained easily in practical applications. Chatzi and Smyth ³⁵ compared the UKF algorithm with the particle filter (PF) method by testing a 3DOF structure with only the first DOF being nonlinear, and it was shown that UKF is the most computationally efficient. Although a discussion was made about the influence of the initial condition intervals for the PF method, no explanation was provided regarding the selection of the initial target parameters in the UKF method.

Chatzi et al. ³⁶ proposed a modified Bouc–Wen (BW) model to capture an oddly shaped hysteresis behavior obtained in an experiment using UKF. Despite not being far from each other, a number of choices for the initial estimation of the BW parameters were guessed, and it was stated that the effects of these initial guesses are negligible. Another experimental implementation of the online model updating using UKF can be found in Song and Dyke, ¹⁶ where a quasi-static cyclic test is performed for offline identification of the BW parameters. The set of identified parameters is later used in the shake table tests as the initial target parameters. In another study, UKF and adaptive constrained UKF methods are compared experimentally; ³⁷ however, the choice of the initial target parameters is not described. The unknown input UKF was also studied on a large-scale concrete shear wall structure. ³⁸ Although the importance of the initial target parameters was discussed for the numerical tests, the initial parameters for the experimental test were chosen based on previous studies.

The importance of the initial guesses for the target parameters is emphasized in most applications of the UKF; yet, to the authors’ knowledge, no nondestructive and practical solutions have been proposed. Moreover, the coupling effect of the initial state covariance matrix and process and measurement noise covariance matrices are not studied comprehensively. This gap exists in spite of the fact that starting the algorithm with arbitrarily chosen parameters may result in the numerical divergence of the process. In this paper, we aim to provide guidelines regarding an optimal and robust system identification with UKF. The joint effect of the initial target parameters and algorithm noise covariance matrices should be investigated in order to be able to diminish the possibility of the algorithm divergence. We demonstrate the impact of the noise parameters on the outcome of the parameter identification on a toy example. Subsequently, we propose a Bayesian optimization method to identify the optimal set of the noise parameters given different sets of initial target parameters. In order to apply this framework to an online parameter identification of a five-story shear frame, a comprehensive numerical analysis of the model is performed. Different sets of initial guesses for the target parameters with different levels of confidence are considered, and their effect on the result of the parameter identification is analyzed. To model the nonlinear behavior of the structure, we use the phenomenological BW model, which has been shown to be a suitable model in this area.^12,39,40 Finally, the application of the proposed framework is validated with an experimental parameter identification of a five-story steel frame.

The remainder of this article is organized as follows: “Uncented Kalman filter in system identification” section briefly describes the UKF formulation and the steps required for a successful implementation; “Optimization of UKF initial conditions” section explains the proposed method for the identification of optimum initial conditions for the UKF algorithm; “Premilinary steps for the implementation of the UKF in shake table tasting” section shows the implementation of the proposed steps to a five-story experimental model; “Parameter identification of the experimental model with UKF” section discusses the results of the UKF implementation on the experimental model with and without taking the preliminary steps; and the final section states the concluding remarks.

Unscented Kalman filter in system identification

Kalman filter is commonly used for state estimation of systems given noisy observations. To identify the hysteretic behavior of a structural system based on a parametric model, the set of target parameters is appended to the state vector. Considering the system state parameters in vector $\mathbf{U}\left(t\right)$ and the structural target parameters in vector $\varphi$ , the continuous state-space formulation for the system identification of a nonlinear dynamic system can be written in the following form

$$\dot{\mathbf{X}}=\mathbf{f}\left(\mathbf{X}\left(t\right), \mathbf{l}\left(t\right),\mathrm{ }\mathbf{w}\left(t\right)\right)$$

(1)

$$\mathbf{Y}\left(t\right)=\mathbf{h}(\mathbf{X}\left(t\right),\mathbf{v}\left(t\right))$$

(2)

where the state-vector $\mathbf{X}\left(t\right)=\{\mathbf{U}\left(t\right), \varphi \}$ , $\mathbf{Y}\left(t\right)$ represents the measurements, l is the input loading, and w and v are system process and measurement noise vectors, respectively, and are assumed to be Gaussian white noise with zero mean. The state equations can be converted to discrete format as follows

$${\mathbf{X}}_k=\mathbf{F}\left({\mathbf{X}}_{k-1},\mathrm{ }{\mathbf{l}}_{k-1}\right)+{\mathbf{w}}_k$$

(3)

$${\mathbf{Y}}_k=\mathbf{H}\left({\mathbf{X}}_k\right)+{\mathbf{v}}_k$$

(4)

In the rest of the article, the covariance matrices for the process and measurement noise vectors are shown by Q and R, respectively. Function F can be obtained from the following equation

$$\mathbf{F}\left({\mathbf{X}}_k,{\mathbf{l}}_k\right)={\mathbf{X}}_k+{\int }_{k\Delta t}^{(k+1)\Delta t}\mathit{f}\left(\mathbf{X}\left(t\right),\mathbf{l}\left(t\right)\right)\mathrm{ }\mathrm{d}t$$

(5)

where $\Delta t$ is the time step. The integration can be solved by numerical methods such as the fourth-order Runge–Kutta method. In each time step, the sigma points are computed for the state parameters to predict the a priori state vector. The sigma points are chosen such that the weighted average and covariance of the propagated points capture the true mean and covariance of the state in the next step. The sigma points can be calculated as

$${\chi }_k=\left[{\widehat{\mathbf{X}}}_k{\widehat{\mathbf{X}}}_k+\sqrt{(L+\lambda ){\mathbf{P}}_k}{\widehat{\mathbf{X}}}_{\mathrm{k}}-\sqrt{(L+\lambda ){\mathbf{P}}_k}\right]$$

(6)

where $L$ is the length of the state-vector X and $\lambda$ is a scaling factor that can be determined as described in Kandepu et al. ⁹ With the calculation of the sigma points, the time update step can be initiated. In this stage, each sigma point is propagated using the nonlinear function F obtained from equation (5). The a priori state estimation and its covariance can then be computed by weighted averaging of the propagated sigma points

$${\left({\mathit{\chi }}_{k+1|k}\right)}_i=\mathbf{F}\left(({\mathit{\chi }}_k{)}_i,{\mathbf{l}}_k\right)$$

(7)

$${\widehat{\mathbf{X}}}_{k+1}^{-}={\sum }_{i=0}^{2L}{W}_i^{\left(m\right)}{\left({\mathit{\chi }}_{k+1|k}\right)}_i$$

(8)

$${\mathbf{P}}_{k+1}^{-}={\sum }_{i=0}^{2L}{W}_i^{\left(c\right)}\left[{\left({\mathit{\chi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{X}}}_{k+1}^{-}\right]{\left[{\left({\mathit{\chi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{X}}}_{k+1}^{-}\right]}^T+\mathbf{Q}$$

(9)

where the dash superscript denotes the a priori stage. Likewise, using the predicted state vectors, the predicted measurements associated with each sigma point and its covariance matrix can be computed

$${\left({\mathit{\psi }}_{k+1|k}\right)}_i=\mathbf{H}\left(({\mathit{\chi }}_{k+1|k}{)}_i,{\mathbf{l}}_{k+1}\right)$$

(10)

$${\widehat{\mathbf{Y}}}_{k+1}^{-}={\sum }_{i=0}^{2L}{W}_i^{\left(m\right)}{\left({\mathit{\psi }}_{k+1|k}\right)}_i$$

(11)

$${\mathbf{P}}_{k+1}^{YY}={\sum }_{i=0}^{2L}{W}_i^{\left(c\right)}\left[{\left({\mathit{\psi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{Y}}}_{k+1}^{-}\right]{\left[{\left({\mathit{\psi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{Y}}}_{k+1}^{-}\right]}^T+\mathbf{R}$$

(12)

where the proper choice of $W_i^{\left(c\right)}$ and $W_i^{\left(m\right)}$ parameters can be found at Kandepu et al. ⁹ Although the noise matrices are added to the state and measurement matrices directly, an augmented state vector can also be considered to propagate the noise through the nonlinear equations as in Kandepu et al. ⁹ The augmented state vector increases the computational costs, and therefore, we preferred to use the former method here.

Next, the measurement update step can be executed using the Kalman gain factor. This factor determines the ratio by which the process calculations and sensor measurements are weighted and applied to the calculation of the next step state vector

$${\mathbf{P}}_{k+1}^{XY}={\sum }_{i=0}^{2L}{W}_i^{\left(c\right)}\left[{\left({\mathit{\chi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{X}}}_{k+1}^{-}\right]{\left[{\left({\mathit{\psi }}_{k+1|k}\right)}_i-{\widehat{\mathbf{Y}}}_{k+1}^{-}\right]}^T$$

(13)

$${\mathit{\zeta }}_{k+1}={\mathbf{P}}_{k+1}^{XY}{\left({\mathbf{P}}_{k+1}^{YY}\right)}^{-1}$$

(14)

$${\widehat{\mathbf{X}}}_{k+1}={\widehat{\mathbf{X}}}_{k+1}^{-}+{\mathit{\zeta }}_{k+1}\left({\mathbf{Y}}_{k+1}-{\widehat{\mathbf{Y}}}_{k+1}^{-}\right)$$

(15)

$${\mathbf{P}}_{k+1}={\mathbf{P}}_{k+1}^{-}-{\mathit{\zeta }}_{k+1}{\mathbf{P}}_{k+1}^{YY}{\mathit{\zeta }}_{k+1}^T$$

(16)

where $\mathit{\zeta }$ is the Kalman gain matrix.

Prior to implementing the UKF process, the initial state vector ${\widehat{\mathbf{X}}}_0$ , the corresponding initial state covariance matrix ${\mathbf{P}}_0$ , the process noise covariance matrix Q, and the measurement noise covariance matrix R should be determined. The selection of suitable entries for these matrices is crucial to the performance of the algorithm. The proposed Bayesian approach for optimal selection of these parameters is described in the next section.

Optimization of UKF initial conditions

During the UKF implementation, the process and measurement noise significantly affect the performance of the parameter identification. ²⁶ These noise values determine the rate of convergence in the state parameters. Diagonal elements of matrix Q provide constant values which are added to the state covariance matrix in each step (equation (9)) to prevent a null Kalman gain. These values are non-zero for state elements that have no constant true values, such as the system parameters in vector $\mathbf{U}$ , and are often considered zero for the target parameters. ³⁰ However, in parameter identification studies, we aim to fit a parametric model to experimental components, and we may not expect a single set of values to be the only correct answer for the defined target parameters. Therefore, for better robustness, we also consider non-zero entries for the diagonal elements of Q that correspond to the elements of $\varphi$ . The process noise vector has the same units as the state vector, and to simplify the formulations, we formed it as

$$\mathbf{w}=\left[\begin{array}{c}\mathbf{1}\\ {\varphi }_0\end{array}\right]$$

(17)

where vector 1 has the same size as $\mathbf{U}$ and ${\varphi }_0$ contain the initial values considered for the target parameters. The process noise covariance matrix can then be tuned by considering a single parameter $w_f$ as

$$\mathbf{Q}=10^{-w_f}\mathbf{I}{\mathbf{w}}^2$$

(18)

where $\mathbf{I}$ is the identity matrix. The noise parameters are assumed to be Gaussian white noise vectors, and thus, the off-diagonal elements of the covariance matrices are zero. Similarly, P₀ is a diagonal matrix where each element reflects the initial uncertainty of the state parameters. ³⁰ Therefore, since the uncertainty for the initial vector $\mathbf{U}$ is zero and we expect the level of uncertainty in the parameters of ${\varphi }_{\hspace{0.17em}}$ to be correlated with the initial guesses for these parameters, the following matrix can be considered for P₀

$${\mathbf{P}}_0=10^{-p_f}\mathbf{I}\left[\begin{array}{c}0\\ {\varphi }_0\end{array}\right]$$

(19)

At the same time, matrix R compensates for the errors in the measurement due to sensor noise and error. Therefore, consistent with the formulations of the Q and ${\mathbf{P}}_0$ , the measurement noise covariance matrix can be formed as

$$\mathbf{R}=\mathbf{I}{\mathbf{v}}^2$$

(20)

where $\mathbf{v}$ is a vector containing the standard deviation of the sensor noise for all measurement channels. Measurement noise is calculated according to the quality of data acquisition. Lower levels of measurement noise result in higher Kalman gains and subsequent increases in the rate of convergence. One should be careful not to choose the measurement noise parameters too small since it opens up the possibility of divergence in the process.

Aside from the noise parameters discussed above, the performance of UKF is affected by the initial guesses about the target parameters, ${\varphi }_0$ . An illustrative example is provided in the next section.

The influence of the initial parameters on the UKF performance

In this section, the influence of the ${\varphi }_0$ parameters is shown by applying UKF to a simplified parameter identification problem. In this example, UKF is utilized to identify the target parameters $a$ and $b$ from the noisy measurements of the function $y=a\sin \left(bt\right)$ . The response is obtained by considering $a=2$ , $b=3$ , and $\mathbf{t}=0:0.01:10$ . To visualize the effect of the initial parameters, a range of different $w_f$ and $p_f$ as well as three initial guesses with different levels of accuracy about the target parameters are considered (Table 1). Gaussian white noise with zero mean and standard deviation of ${10}^{-2}$ is added to the observation signal in all cases.

Table 1.

Different initial guesses for the illustrative example.

	Case 1: true	Case 2: close	Case 3: inaccurate
${\varphi }_0$ parameters	$a=2$ , $b=3$	$a=2.4$ , $b=2.4$	$a=0.2$ , $b=15$

The results of the parameter identification are evaluated by calculating the root mean square (RMS) error of the distance from the identified parameters to the true values of the target parameters. Figure 1 shows the RMS values for different noise and initial parameters. It can be observed that, in all three cases, we only obtain accurate results for certain combinations of noise parameters. If $w_f$ and $p_f$ are chosen too large, the UKF process practically ends up delivering the same parameters as the initial guesses, and hence we observe the accurate predictions at the top right corner of the Case 1. Moreover, as the initial guesses for the state parameters draw apart from the true values in Cases 2 and 3, the region in which we obtain accurate results shrinks. This observation emphasizes the need for carefully selecting the initial noise parameters in situations where likely guesses about the initial target parameters are difficult to make.

Figure 1.

Results of the UKF implementation on the illustrative example with different initial conditions. (a) Case 1, (b) Case 2, and (c) Case 3.

The proposed optimal initialization process for UKF

Based on the observations from the illustrative example in the previous section, it can be inferred that the need for an inclusive approach determine the initial conditions of the UKF algorithm is crucial. Therefore, in this manuscript, a framework that is designed to ensure a successful UKF parameter identification is proposed.

Parameter identification is usually performed with phenomenological models where parameters do not convey physical meanings. Therefore, the proposed method begins with creating a simple constitutive model, such as a bilinear model, to get a good sense regarding the initial guesses about the target parameters. This model is created with the knowledge about the basic properties of the experimental specimen including geometric dimensions, yield stress, mass, etc. Then, a phenomenological model is fitted to the numerical constitutive model in an offline process. The identified parameters for this model are stored to be used as the initial guesses for the final experimental testing. Next, the virtual emulation of the experiment is performed by considering the fitted phenomenological model as the actual experimental specimen and intentionally distancing the guesses about the initial target parameters from the identified parameters. These inaccurate guesses are planned to account for the uncertainty of the guessed initial target parameters in the actual experiment. By emulating the experiment, the noise combinations that perform well for inaccurate initial target parameters can be identified. The sets of combinations can be divided into different zones. Based on the confidence in the proximity of the fitted model to the experimental specimen, the operator can select either a set of noise parameters from an “Accurate” zone or a “Robust” zone. Parameters in the Accurate zone are intended to deliver the best UKF performance, provided that the initial guesses are chosen close to the actual parameters. On the contrary, the Robust zone provides noise parameters that prevent the divergence of the algorithm for incautiously guessed initial target parameters. Since every condition in the experimental test is considered in the simulation, the authors hypothesize that the noise combinations found for the emulated model will adequately fit the experiment.

The identification of the suitable noise parameters requires tedious virtual implementations of the UKF process on the initially fitted phenomenological model. Therefore, to reduce the computational costs, a Bayesian optimization algorithm is utilized to find the mentioned zones. Details of the Bayesian optimization are described in the next section.

Bayesian optimization of UKF algorithm parameter

It has been customary to obtain a reasonable set of noise parameters $\mathit{\theta }=[\mathbf{Q}, \mathbf{R}, {\mathbf{P}}_0]$ by trial and error in the current literature. To address the existing gap, the authors propose a Bayesian optimization with the Gaussian process approximation^41,42 strategy to find the best set of parameters in this study.

We formulate Bayesian optimization to find the global minimum of the objective function over a compact set C , which consists of a range of different combinations for the parameters in $\mathit{\theta }.$ The RMS value of the difference between the measured displacement and the displacement obtained from the final identified BW parameters is considered as the objective function. This function can be optimized with the following equation

$$\underset{\mathit{\theta }\in \mathit{C}}{\text{min}\mathrm{ }}g\text{(}\mathit{\theta }\text{)}$$

(21)

where $g$ calculates the RMS value. Since the problem is non-parametric, the objective function should be queried at all instances of Q, R, P₀ in order to find the extremum. However, each query is equivalent to performing the computationally expensive UKF emulation. In contrary to the grid search approach, Bayesian optimization provides the best “next sample” in order to reduce the number of total function evaluations.

This procedure starts with choosing a few instances randomly and fitting Gaussian processes to the observations. Then, an exploration–exploitation trade-off should be considered for the choice of the next sampling point. The exploration contributes to areas with high uncertainty, and the exploitation takes advantage of the points that are more likely to be closer to the extremum. The trade-off can be managed in the form of an improvement-based acquisition function. The Bayesian approach suggests using the following acquisition function which automatically balances the trade-off

$$EI\left(\mathit{\theta }\right)=\left\{\begin{array}{l}\left(\mu \left(\mathit{\theta }\right)-g\left({\mathit{\theta }}^{+}\right)\right)\mathrm{\Omega }\left(Z\right)+\sigma \left(\mathit{\theta }\right)\omega \left(Z\right)\text{ if }\sigma \left(\mathit{\theta }\right)>0\\ 0\text{ if }\sigma \left(\mathit{\theta }\right)=0\end{array}\right.$$

(22)

$$Z=\frac{\mu \left(\mathit{\theta }\right)-g\left({\mathit{\theta }}^{+}\right)}{\sigma \left(\mathit{\theta }\right)}$$

(23)

where $\mu$ is the mean and $\sigma$ is the standard deviation of the Gaussian processes at $\theta$ . Operators $\mathrm{\Omega }$ and $\omega$ represent the cumulative distribution function and PDF of the parameters, respectively, and ${\mathit{\theta }}^{+}={\mathrm{argmax}}_{{\mathit{\theta }}_i\in {\mathit{\theta }}_{1:t}}g\left({\mathit{\theta }}_i\right)$ . The best “next sample” point can be obtained by maximizing the EI function. More details and derivation of equations can be found at Brochu et al. ⁴³

For the sake of simplicity and better visualization of the results, the measurement noise parameters are considered constant in this article and are determined based on the measured noise from the sensors. However, one can include $\mathbf{R}$ in the Bayesian optimization for higher accuracy. The proposed steps for achieving a set of optimum preliminary conditions for the final experiment are summarized in the flowchart shown in Figure 2.

Figure 2.

Flowchart of the proposed approach.

Preliminary steps for the implementation of the UKF in shake table testing

The proposed approach is implemented in the experimental parameter identification of a five-story shear frame. The preliminary steps prior to the experimental testing are described in this section, and the experimental procedure is discussed in “Parameter identification of the experimental model with UKF” section. In the following subsections, the state-space formulation for parameter identification of an MDOF BW model is derived, and the required steps before the final experiment follow next.

Equation of motion for a nonlinear MDOF system

As stated in the “Introduction” section, in this study, the nonlinear behavior of steel material is modeled by a classical BW model, and degradation and pinching effects are neglected. The equation of motion for an MDOF nonlinear system can be considered as follows

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{S}=\mathbf{l}\left(t\right)$$

(24)

where M and C are the mass and damping matrices, and u, $\dot{\mathbf{u}}$ , and $\ddot{\mathbf{u}}$ are displacement, velocity, and acceleration vectors, respectively. For the classical BW model, the restoring force of the ith story S_i can be formed as

$$S_i\left(\mathbf{u},\mathbf{z}\right)={\alpha }_i{k}_i(u_i-u_{i-1})+\left(1-{\alpha }_i\right)k_i{z}_i$$

(25)

where

$${\dot{z}}_i=A_i({\dot{u}}_i-{\dot{u}}_{i-1})-{\beta }_i\left|{\dot{|u}}_i-{\dot{u}}_{i-1}\right|||{\left|z_i|\right|}^{n-1}{z}_i-{\gamma }_i({\dot{u}}_i-{\dot{u}}_{i-1})|{\left|z_i|\right|}^n$$

(26)

where $z_i$ is the hysteresis component, $k_i$ is the stiffness coefficient, and ${\alpha }_i\in [0, 1]$ indicates the level of the nonlinearity of the system. A_i, ${\beta }_i$ , n, and ${\gamma }_i$ are the coefficients determining the shape of the hysteresis loops. The effect of each parameter on the shape of the hysteresis cycles can be found in Chatzi et al. ³⁶ and Ou ⁴⁴

Considering the equation of motion for a q-story shear building with identical lumped masses m located at floors, the mass matrix can be considered as $m{\mathbf{I}}_{q\times q}$ and the damping matrix as

$${\left[\begin{array}{ccccc}{c}_1+c_2& -c_2& \cdots & \mathrm{ }& 0\\ \mathrm{ }& \ddots & \mathrm{ }& \mathrm{ }& \mathrm{ }\\ ⋮& -c_i& c_i+c_{i+1}& -c_{i+1}& ⋮\\ \mathrm{ }& \mathrm{ }& \mathrm{ }& \ddots & \mathrm{ }\\ 0& \mathrm{ }& \cdots & -c_q& c_q\end{array}\right]}_{q\times q}$$

(27)

Assuming identical components, the damping matrix can be formed using only two parameters $c_1$ and $c_2$ as

$${\left[\begin{array}{ccccc}{c}_1+c_2& -c_2& \cdots & \mathrm{ }& 0\\ \mathrm{ }& \ddots & \mathrm{ }& \mathrm{ }& \mathrm{ }\\ ⋮& -c_1& c_1+c_2& -c_1& ⋮\\ \mathrm{ }& \mathrm{ }& \mathrm{ }& \ddots & \mathrm{ }\\ 0& \mathrm{ }& \cdots & -c_1& c_2\end{array}\right]}_{q\times q}$$

(28)

Therefore, the equations of motion can be decoupled as follows

$$m\ddot{u_1}+(c_1+c_2){\dot{u}}_1-c_2{\dot{u}}_2+S_1-S_2=l_1\left(t\right)$$

$$⋮$$

$$m\ddot{u_i}+{(c}_1+c_2){\dot{u}}_i-c_1\left({\dot{u}}_{i-1}+{\dot{u}}_{i+1}\right)+S_i-S_{i+1}=l_i\left(t\right)$$

(29)

$$⋮$$

$$m\ddot{u_q}+c_2{\dot{u}}_q-c_1{\dot{u}}_{q-1}+S_q=l_q\left(t\right)$$

Equation set (29) can be used for the finite element model updating of a nonlinear shear building.

State space formulation of the MDOF system

In order to formulate the dynamic response problem of an MDOF structure into nonlinear parameter identification formulation, the state-space form of the MDOF equation of motion is constructed. The floor accelerations are considered as measurements, while the ground acceleration signal is utilized as the input control to the system. The state vector consists of the system-state variables followed by the BW parameters

$$\mathbf{X}=\left\{\mathbf{U}\left(t\right),\mathrm{ }\varphi \right\}=\left\{\left\{{\mathbf{u}}_{1\times q},{\dot{\mathbf{u}}}_{1\times q},{\mathbf{z}}_{1\times q}\right\},\left\{{\mathbf{A}}_{1\times q},{\mathbf{\beta }}_{1\times q},{\mathbf{\gamma }}_{1\times q},{\mathbf{\alpha }}_{1\times q},{\mathbf{c}}_{1\times q},{\mathbf{k}}_{1\times q}\right\}\right\}$$

(30)

The nonlinear state-space equation can be written as

$$\dot{\mathbf{X}}=\mathbf{f}\left(\mathbf{X}\left(t\right),\mathbf{l}\left(t\right)\right)$$

(31)

In the case of seismic loading, the input $\mathbf{l}$ is defined as

$$\mathbf{l}\left(t\right)=-{\ddot{u}}_{g}m\mathbf{r}$$

(32)

where ${\ddot{u}}_g$ is the ground acceleration and r is the displacement transformation vector. Considering equations (29), equation (31) can be written as

$$\mathbf{f}\left(\mathbf{X}\left(t\right),\mathbf{l}\left(t\right)\right)=\left[\begin{array}{c}\dot{\mathbf{u}}\\ \ddot{\mathbf{u}}\\ \begin{array}{c}\dot{\mathbf{z}}\\ \dot{\varphi }\end{array}\end{array}\right]=\mathrm{ }\left[\begin{array}{c}\dot{\mathbf{u}}\\ m^{-1}\left(l_1\left(t\right)-\left(c_1+c_2\right){\dot{u}}_1+c_2{\dot{u}}_2-k_1\left({\alpha }_1{u}_1+\left(1-{\alpha }_1\right)z_1\right)+k_2\left({\alpha }_2(u_2-u_1)+(1-{\alpha }_2)z_2\right)\right)\\ \begin{array}{c}⋮\\ m^{-1}\left(l_i\left(t\right)-c_i\left({\dot{u}}_i{-\dot{u}}_{i-1}\right)-c_{i+1}\left({\dot{u}}_i{-\dot{u}}_{i+1}\right)-k_i\left({\alpha }_i{(u}_i-u_{i-1})+\left(1-{\alpha }_i\right)z_i\right)+k_{i+1}\left({\alpha }_{i+1}(u_{i+1}-u_i)+(1-{\alpha }_{i+1})z_{i+1}\right)\right)\\ \begin{array}{c}⋮\\ m^{-1}\left(l_q\left(t\right)-c_q\left({\dot{u}}_q{-\dot{u}}_{q-1}\right)-k_q\left({\alpha }_q{(u}_q-u_{q-1})+\left(1-{\alpha }_q\right)z_q\right)\right)\\ \begin{array}{c}{A}_1{\dot{u}}_1-{\beta }_1\left|\hspace{0.17em}|{\dot{u}}_1|\hspace{0.17em}\right|{\left|\hspace{0.17em}|z_1|\hspace{0.17em}\right|}^{n-1}{z}_1-{\gamma }_1{\dot{u}}_1{\left|\hspace{0.17em}|z_1|\hspace{0.17em}\right|}^n\\ ⋮\\ \begin{array}{c}{A}_i({\dot{u}}_i-{\dot{u}}_{i-1})-{\beta }_i\left|\hspace{0.17em}|{\dot{u}}_i-{\dot{u}}_{i-1}|\hspace{0.17em}\right|{\left|\hspace{0.17em}|z_i|\hspace{0.17em}\right|}^{n-1}{z}_i-{\gamma }_i({\dot{u}}_i-{\dot{u}}_{i-1}){\left|\hspace{0.17em}|z_i|\hspace{0.17em}\right|}^n\\ ⋮\\ \begin{array}{c}{A}_q({\dot{u}}_q-{\dot{u}}_{q-1})-{\beta }_q\left|\hspace{0.17em}|{\dot{u}}_q-{\dot{u}}_{q-1}|\hspace{0.17em}\right|{\left|\hspace{0.17em}|z_q|\hspace{0.17em}\right|}^{n-1}{z}_q-{\gamma }_q({\dot{u}}_q-{\dot{u}}_{q-1}){\left|\hspace{0.17em}|z_q|\hspace{0.17em}\right|}^n\\ 0_{6q\times 1}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$$

(33)

At the same time, the acceleration measurement vector can be calculated as follows:

$$\mathbf{Y}=\mathbf{h}\left(\mathbf{X}\left(t\right),\mathbf{l}\left(t\right)\right)$$

(34)

where

$$\mathbf{Y}=\ddot{\mathbf{u}}-\frac{\mathbf{l}\left(t\right)}{m}=\left[\begin{array}{c}{m}^{-1}\left(-\left(c_1+c_2\right){\dot{u}}_1+c_2{\dot{u}}_2-k_1\left({\alpha }_1{u}_1+\left(1-{\alpha }_1\right)z_1\right)+k_2\left({\alpha }_2(u_2-u_1)+(1-{\alpha }_2)z_2\right)\right)\\ ⋮\\ m^{-1}\left(-c_i\left({\dot{u}}_i{-\dot{u}}_{i-1}\right)-c_{i+1}\left({\dot{u}}_i{-\dot{u}}_{i+1}\right)-k_i\left({\alpha }_i{(u}_i-u_{i-1})+\left(1-{\alpha }_i\right)z_i\right)+k_{i+1}\left({\alpha }_{i+1}(u_{i+1}-u_i)+(1-{\alpha }_{i+1})z_{i+1}\right)\right)\\ ⋮\\ m^{-1}\left(-c_q\left({\dot{u}}_q{-\dot{u}}_{q-1}\right)-k_q\left({\alpha }_q{(u}_q-u_{q-1})+\left(1-{\alpha }_q\right)z_q\right)\right)\end{array}\right]$$

(35)

In case the same components are used for all stories of the shear model, the vectors of BW parameters in vector $\varphi$ can be substituted with scaler values.

Experimental model

The experimental model represents a five-story shear building with steel columns. This specimen is manufactured with identical floor masses and column components. The details of the floor sections and columns are shown in Figure 3. The contributing mass of each story is 23.2 kg, and the total linear stiffness of the components in each floor is approximated as 1.82e5 kN/m.

Figure 3.

The five-story building model.

At this stage, the measurement noise levels are determined in idling situations. The standard deviations of the five acceleration sensors were measured in units of $m/s^2$ as $\mathbf{v}=\{0.0024, 0.0023, 0.0018, 0.0013, 0.0022\}$ .

Identification of the ${\varphi }_0$ parameters

A preliminary numerical analysis is required to obtain reasonable guesses for the ${\varphi }_0$ parameters and optimal set of noise parameters $\mathbf{Q}$ and ${\mathbf{P}}_0$ for the experimental test. Therefore, a bilinear model with kinematic hardening is developed with regards to the geometry and materials of the experimental setup. Uniaxial cold and hot tension tests were performed on the steel material batch used in the lab manufacturing. The numerical model uses a simplified bilinear Steel01 material, which is created in OpenSees with a yield stress of 460 MPa and the modulus of elasticity of 200 GPa. The strain hardening ratio is assumed to be 5% of the initial elastic tangent. The enforced displacement time history and the cyclic behavior of the numerical model can be seen in Figure 4.

Figure 4.

The behavior of the constitutive bilinear model: (a) displacement input and (b) hysteretic cycles.

Based on the numerical response of the 1-DOF system, a set of BW parameters are fitted to the bilinear model, as shown in Figure 5(a). In order to simulate the uncertainty of initial state parameters in the experimental test, inaccurate ${\varphi }_0$ parameter sets are also considered. The two other ${\varphi }_0$ vectors are tagged as “Confident estimation” (Figure 5(b)) and “Inaccurate estimation” (Figure 5(c)). The second case represents a rough estimation of the actual parameters, which might be practical for complex models. And the third case represents a careless selection of the parameters that create reasonable hysteresis cycles despite being far from the actual case. The full set of the considered ${\varphi }_0$ vectors for each case is shown in Table 2.

Figure 5.

The comparison between cyclic behaviors of the considered BW models with the bilinear model. (a) Case 1, (b) Case 2, and (c) Case 3.

Table 2.

The values of the initial guesses for the target parameters of each case.

${\varphi }_0$ parameters	Case 1: true	Case 2: confident	Case 3: inaccurate
A	0.84	0.67	1.85
$\mathit{\beta }$	1395.21	1116.2	2790.41
$\mathit{\gamma }$	40.44	32.35	55.00
$\mathit{\alpha }$	0.09	0.06	0.09
${\mathit{c}}_1$	73.68	58.94	221.04
${\mathit{c}}_2$	91.11	72.88	27.33
K	182,000.00	145,600.00	18,200.00

The parameter n in a BW model is an exponential constant, and a minor variance may result in the divergence of the algorithm. This issue worsens in the presence of signal noise, and it is beneficial to consider $n$ as a constant in practice. ⁴⁵ Therefore, obtained from the offline fitting of the BW parameters to the numerical bilinear model, $n$ is considered to be 1.65 and will not be considered a target parameter.

Optimization of the UKF noise parameters

To identify a set of optimal noise parameters for the experimental test, different combinations are tested on emulated test setups, and different zones are identified with Bayesian optimization. First, a time history analysis with the specified ground motion is performed on the fitted BW model (Case 1 in Table 2), and acceleration responses are recorded. Then, different rounds of UKF are performed using different ${\varphi }_0$ parameters of Table 2 as initial guesses for the target parameters and the acceleration responses of the Case 1 as measurements. The excitation is 30 s of a scaled El-Centro earthquake both in time and amplitude with a maximum amplitude of 0.65 g and with a time step of 0.005 s (Figure 6).

Figure 6.

The time series of the El-Centro earthquake excitation.

The measured signals in the virtual emulation are contaminated with white Gaussian noise at the same level as the measured noise from actual sensors in vector $\mathbf{v}$ . The matrix R is then calculated according to equation (20). In order to find the best noise combinations for matrices $\mathbf{Q}$ and ${\mathbf{P}}_0$ , optimal sets of $w_f$ and $p_f$ are solved with a two-dimensional Bayesian optimization, as explained in “Bayesian optimization of UKF algorithm parameter” section. The considered domain for $w_f$ ranges from 0 to 25 with 0.5 increments, and $p_f$ is chosen to range from –1 to 10, totaling 600 combinations. To start the Bayesian optimization process, 30 combinations are randomly chosen, and the performance of the parameter identification is evaluated. Then, these points are considered as the input to the Gaussian processes, and the next sampling point is calculated from equation (22). The process continues until a convergence to the minimum value achieves.

The outcome of each optimization process is presented in a 2D domain showing all the combinations of the noise parameters considered. In each figure, the yellow-colored spots show the initial 30 points. The other spots marked as green indicate the points suggested by Bayesian optimization for evaluation, and the surface shows the mean value for the fitted Gaussian process model. For these surfaces, areas with lower RMS errors are illustrated with darker colors. Spots marked with a red symbol denote the divergence of the process in that noise combination. The optimization solver is trained to avoid the red spots in order to focus on areas where the algorithm converges with better performance.

The optimization path is shown in terms of RMS value in Figure 7 and also in Figure 8 as a 2D map. Figure 7 shows the suggestions of the Bayesian optimization through the convergence to the baseline minimum, which is the global optimum in the chosen domain. Comparing the minimum identification errors (the red-dashed lines in Figure 7) among the three cases, it can be inferred that the minimum error found significantly increases with a worse choice of the initial target parameters.

Figure 7.

Illustration of individual sample points, and the path suggested by Bayesian optimization. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 8.

The domain considered for $w_f$ and $p_f$ and the optimization results for different cases. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 8 shows that a Robust zone can be described by a triangular shape for all three cases. Although the performance is better for Cases 1 and 2, in this zone, the process converges regardless of the initial target parameters. This zone is enclosed with cyan dash-dotted lines in Figure 8. Moreover, a smaller trapezoid area can be marked for the first and second cases where the error is significantly lower than other areas and can be considered as the Accurate zone. It should be pointed out that choosing a noise combination in this area will result in the divergence of the process for Case 3 where the initial target parameters are chosen carelessly. The Accurate zone is enclosed with magenta-colored dashed lines in Figure 8.

Parameter identification of the experimental model with UKF

The experimental model is mounted on a 6-DOF shake table in the Intelligent Infrastructure System Lab (IISL; https://engineering.purdue.edu/IISL/) at Purdue University for the dynamic testing. In this setup, the El-Centro ground motion is imposed laterally to the building’s axis. Accelerometers and optical sensors are used to measure the absolute responses of the frame, as indicated in Figure 9. The data acquisition system recorded the structural acceleration responses with a sampling frequency of 2048 Hz. Also, the displacement responses are measured using a 6D krypton optical tracking system, capturing the position of LED sensors on model’s base and floors. The krypton system recorded the responses with a sampling frequency of 60 Hz.

Figure 9.

Shake table test configuration.

To test the UKF performance with optimized algorithm parameters determined from the numerical study, we only need to use the target parameters of Case 1 in Table 2 as ${\varphi }_0$ and a set of noise parameters from the identified zones. However, to show the effect of the preliminary steps on the outcome of the parameter identification of the experimental model, the UKF process is also performed with different initial conditions. First, a noise combination is chosen from the Accurate zone for the first two cases of ${\varphi }_0$ parameters with $w_f=14$ and $p_f=7$ , and second, a combination from the Robust zone for the third case of the initial state parameters as $w_f=8$ and $p_f=8$ . Displacement responses of the experimental tests are shown in the form of the time history responses of the first and fifth stories of the model. Moreover, the responses of the numerical simulations with initial state parameters and the parameters that the UKF process converged to are plotted over the experimental signals. The displacement signals from the initial guesses and identified parameters are compared to the experimental signal for Case 1 in Figures 10 and 11, Case 2 in Figures 12 and 13, and Case 3 in Figures 14 and 15. Also, the response spectrums in frequency domain and the force-displacement graphs of the first story components for all cases are shown in Figures 16 and 17, respectively.

Figure 10.

Displacement response of the first story of the model for Case 1 of the initial BW parameters.

Figure 11.

Displacement response of the fifth story of the model for Case 1 of the initial BW parameters.

Figure 12.

Displacement response of the first story of the model for Case 2 of the initial BW parameters.

Figure 13.

Displacement response of the fifth story of the model for Case 2 of the initial BW parameters.

Figure 14.

Displacement response of the first story of the model for Case 3 of the initial BW parameters.

Figure 15.

Displacement response of the fifth story of the model for Case 3 of the initial BW parameters.

Figure 16.

The amplitude frequency spectrums of the response displacement. The top row shows the responses of the first story, and the bottom row shows the fifth story responses. (a) Case 1, (b) Case 2, and (c) Case 3.

Figure 17.

The hysteresis behavior of the first story component under seismic loading. (a) Case 1, (b) Case 2, and (c) Case 3.

The quantitative improvements obtained from UKF identified parameters over the initial ${\varphi }_0$ parameters are shown in Table 3. Three different criteria are defined to measure the performance, namely, the difference between the time histories and power spectrums in terms of the RMS error, and the dominant frequency of the displacement response. The values in Table 3 are averaged over all stories of the model. Some values are not available for Case 3 because of the divergence of the time history analysis for the initial parameter set. From Figures 10 and 11, it can be inferred that the initial parameter set for Case 1 can accurately imitate the responses of the actual experimental test. Therefore, we can infer that fitting a BW model to the assumed constitutive model has guided us to a very close estimation of the best fit of a BW set. Nevertheless, the UKF process has improved the accuracy further according to Table 3, and this improvement can be seen in the magnified plots of Figures 10 and 11. On the other hand, Figures 12 and 13 show that the small disturbance in the initial target parameters of Case 2 can change the dynamic behavior of the model significantly. Similarly, the UKF algorithm converged to a set of BW parameters that follow the experimental signals carefully. Finally, although the time history analysis diverges at the very early stages of the process with ${\varphi }_0$ parameters of Case 3, the UKF process still converges to a set of BW parameters, that, except the damping ratio, can model the behavior of the experimental model (Figures 14 and 15). In fact, same as Cases 1 and 2, the identified dominant frequency for Case 3 is still 98% accurate, according to Table 3. This observation indicates that choosing a noise combination in the Robust zone will result in the convergence of the UKF process, though the eventual results may not be optimal.

Table 3.

Rate of improvement in parameter identification with UKF in different cases.

Case	RMS time history			RMS response spectrum			Signal’s dominant frequency (Hz)
	Initial	Final	improvement	Initial	Final	improvement	Experimental	Initial	Final
1	0.0025	0.0021	16.0%	0.2197	0.1107	49.6%	3.60	3.73	3.67
2	0.0087	0.0021	75.9%	0.6852	0.1101	83.9%	3.60	2.93	3.67
3	∞	0.0039	∞	∞	0.2547	∞	3.60	n/a	3.67

To find out the performance of the UKF algorithm on the experimental testing given other noise combinations, the results of a grid search for all the combinations of the $w_f$ and $p_f$ parameters are shown in Figure 18. The heatmap surface shows the RMS value for each combination. Also, the borders of the zones identified in the numerical analysis and the chosen combinations based on those zones are shown in this figure. It can be seen that the heatmaps show good correlations with the result of the Bayesian optimization in Figure 8. Therefore, the zones are valid in the experimental analysis as well, and the chosen combinations for the experimental testing performed earlier are well within the aimed areas.

Figure 18.

Performance of the UKF algorithm on the experimental test setup given different initial target parameters. The borders of the robust and accurate zones, as well as the chosen combinations for the experimental testing are shown on top of the experimental results. (a) Case 1, (b) Case 2, and (c) Case 3.

Finally, it can be argued that the Bayesian optimization of the UKF parameters can help preventing the divergence of the process in the expensive experimental tests. Also, if good chances of convergence with a high performance are required, one should obtain a decent estimation of the initial target parameters.

Conclusion

Parameter identification using the UKF algorithm is a tricky task and needs careful attention to the details. Several parameters contribute to the outcome of the algorithm, and in case of experimental testing, a preliminary study should be performed in order to obtain a reasonable set of parameters to start the process with. In this paper, a Bayesian optimization approach for finding the initial parameters of UKF algorithm is proposed and studied. Utilizing the optimized results, a combination of the noise parameters is chosen, and the parameter identification process is implemented on a five-story model, experimentally. Specific findings and contributions are summarized as follows:

The hypothesis is proved that the process noise covariance matrix as well as the initial state covariance matrix induces large variations in the outcome of the UKF algorithm, and poor judgment may result in the divergence of the process.

Given different initial BW model parameters, it is validated that the Bayesian optimization can automatically select the best noise combination required to run the UKF. The findings yield the same optimal identification performance as compared to the combinations determined from a grid search approach.

The optimization approach demonstrates that a trade-off exists between the robustness of the UKF algorithm and the final identification accuracy. Therefore, a suitable choice of the initial state parameters can avoid the divergence of the UKF process, but furthermore, if a high accuracy outcome is desired, choosing a splendid combination for the noise parameters is undeniable.

Meanwhile, the minimum identification error, which is the theoretical best identification performance that can be achieved, is affected by the initial guesses for the target parameters of the BW model. If the initial guesses are closer to the true parameters, the possible minimum error is smaller. Therefore, the selection of initial BW parameters is more important and should be carefully managed.