Identification of structural parameters and boundary conditions using a minimum number of measurement points

Abstract

This article proposes a novel methodology that uses mathematical and numerical models of a structure to build a data set and determine crucial nodes that possess the highest sensitivity. Regression surfaces between the structural parameters and structural output features, represented by the natural frequencies of the structure and local transmissibility, are built using the numerical data set. A description of a possible experimental application is provided, where sensors are mounted at crucial nodes, and the natural frequencies and local transmissibility at each natural frequency are determined from the power spectral density and the power spectral density ratios of the sensor responses, respectively. An inverse iterative process is then applied to identify the structural parameters by matching the experimental features with the available parameters in the myriad numerical data set. Three examples are presented to demonstrate the feasibility and efficacy of the proposed methodology. The results reveal that the method was able to accurately identify the boundary coefficients and physical parameters of the Euler-Bernoulli beam as well as a highway bridge model with elastic foundations using only two measurement points. It is expected that the proposed method will have practical applications in the identification and analysis of restored structural systems with unknown parameters and boundary coefficients.

Keywords

structural model validation / eigenvalue problem / response surface / inverse problems

Introduction

Structural modeling plays an essential role in any model updating or structural health monitoring process [¹]. The structural model is normally initiated based on the nominal values of material properties provided by the manufacturing codes and design drawings as well as the analyst’s assumptions about the connectivity and element types and sizes [²]. It is inevitable that experimental and numerical models will deviate from each other for reasons such as modeling idealizations, model simplifications, and roundoff errors in numerical models associated with computer-based processes [³]. Structural model updating or structural model validation (SMV) can be categorized as the first essential step in practical structural health monitoring [⁴]. During the updating process, the numerical model is tuned to match the in situ real-world model normally by minimizing the discrepancies between their responses.

The most important system parameters in an SMV are the material properties, model geometry, and boundary conditions (BCs). Although there has been extensive focus on updating the system parameters, the updating of BCs has not received similar attention [⁵,⁶]. In most applications, the BCs are idealized as either fixed or free. However, BCs can deteriorate by aging and rust and thus can be a major source of system uncertainty in structural health monitoring and updating processes [⁷,⁸]. While there are guides for the deterioration of material properties, such as Young’s modulus of elasticity and density, due to environmental effects and manufacturing deficiency, there are no such guides for BCs [⁹,¹⁰]. This process can become more complicated when dealing with detailed real-world engineering structures whose optimization schemes require accurate local and global information [¹¹].

The determination of system features is a crucial step in any SMV process. The conventional system features for SMV are the results of an eigenvalue problem (EVP) represented by the natural frequencies (NFs) and mode shapes [¹²]. In a linear system, the relationship between the force and displacement in the static form is linear, whereas the relationship between the system parameters and its features is nonlinear. This could be the main reason that most optimum points returned by customary SMV algorithms are not reliable and are usually local optimum points [¹³,¹⁴]. In addition, there is no way to check the uniqueness of the results to justify the entire process.

The solution of the inverse EVP is one of the most topical subjects in the structural engineering field; it is essential for the SMV process and has been investigated by many authors. Nanthakumar et al. [¹⁵] suggested an algorithm to solve the inverse problem in order to detect material interfaces iteratively using the extended finite element model (FEM). Another SMV method was proposed by Mao and Dai [¹⁶], who used a quadratic inverse eigenvalue solver with incomplete modal data. An inverse structural modification process was developed by Tsai et al. [¹⁷], who solved the inverse problem of receptance instead of the original matrices or modal model. The method iteratively minimizes the dynamic stiffness modification matrix defined in their study and the corresponding objective functions. Anitescu et al. [¹⁸] proposed a method to solve the inverse problem of a complex-valued Helmholtz equation by using artificial neural networks and a grid point generation engine to approximate a differential equation in its domain and the BCs. Experimental SMV studies were performed by Nehete et al. [¹⁹] to obtain accurate information about the vibroacoustic cavity by using an inverse eigen sensitivity method; the authors were able to identify the BC properties of their 3D model.

A major challenge with SMV implementation is determining the most significant structural responses with the smallest number of sensors. In addition, the locations and distributions of the sensors must be appropriately assigned to collect structural features that have sufficient information to identify system parameters [²⁰]. Conventional optimal sensor installation is normally performed based on information entropy and probabilistic tools. The information entropy, as a scalar value, measures the quality of structural information using various clusters of sensor configurations [²¹]. Stochastic simulations using the Bayesian framework must be performed repeatedly to evaluate the optimal sensor locations and uncertainties. An extensive review of the prevalent optimal sensor placement was performed by Mallardo and Aliabadi [²²]. The concurrent methods depend mostly on probabilistic, heuristic, artificial, and genetic algorithm theories [²³,²⁴].

The current structural model updating algorithms have the following disadvantages: (i) system responses normally require numerous degrees of freedoms (DOFs) for mode shape extraction; (ii) the BCs of the system are traditionally assumed to be known and are not considered in the updating process; (iii) there is no way to determine whether the returned value is genuine; (iv) gradient-based and heuristic methods are typically executed with sequential, implicit, and black-box steps that will turn out to be the optimum point closest to the initial guess of the analyst. Several apt treatments are offered to eliminate existent flaws during the process.

A novel approach is presented in this work to identify the uncertain parameters of a system as well as its BCs using the minimum number of measurement points. A new, exquisite method for assigning appropriate weights to the objective function is presented. In this proposed method, a returned point’s stability and uniqueness can be investigated by a visual inspection toolbox, which is unprecedented for SMV analysis. The parallel coordinate toolbox is also presented to enhance the SMV’s confidence in determining whether the returned points are in a unique region. Furthermore, the economic cost of the experiment is reduced by recognizing and then mounting a sufficient number of sensors, while the required number of structural modes are excited.

This article is organized as follows. The methodology section introduces a novel comprehensive algorithm for authentic SMVs. The algorithm uses a multi-stage process to determine a system’s properties with high confidence. Three examples are presented in the example and verification section, ranging from a simple model to more sophisticated ones, to show the feasibility of the proposed method in various situations. The possible experimental implementation section elaborates on the implementation of the proposed method in real-life scenarios. The article ends with the discussion and conclusion section.

Methodology

In this section, the steps required for conventional SMV are presented first. This is followed by a description of the proposed algorithm.

Conventional structural model updating process

In this process, a common objective function, which defines the difference between the analytical and experimental models, is minimized subject to some constraints. The following is an example [¹³,²⁵]:

(1)\min :\hspace{0.17em}J(\underset{\_}{X})={\displaystyle \sum _{i=1}^N{\omega }_i^f}{\left(\frac{f_{a,i}-f_{e,i}}{f_{e,i}}\right)}^2+{\displaystyle \sum _{i=1}^N{\omega }_i^{\varphi }}{\left(\frac{1-MA{C}_i}{MA{C}_i}\right)}^2,

(2)S.T.:\hspace{0.17em}\hspace{0.17em}LB\le \underset{\_}{X}\le UB

(3){\displaystyle \sum _{i=1}^N\left({\omega }_i^f+{\omega }_i^{\varphi }\right)=1},

where \underset{\_}{X}=\{x_{\mathrm{1}}\dots x_M\}is a vector that contains uncertain parameters defined as random variables (RVs). N and M are the total number of modes and total number of uncertain system parameters, respectively. LB and UB are the lower and upper bounds of the RVs, respectively. The feasible region (FR) of the unknown parameters is defined by those bounds. MAC is the modal assurance criterion, which measures the correlation between a pair of local mode shapes [³,²⁶]. Finally, f_{a,i}, f_{e,i}, {\omega }_i^f, and {\omega }_i^{\varphi } are the ith analytical NF, experimental NF, NF weighting coefficient, and mode shape coefficient of the ith mode, respectively.

This formulation normally uses a single objective function that comprises weighted terms, where each weight indicates the relative degree of importance of that term. This process faces a critical pitfall in assigning appropriate weights to each term in order to determine the optimum point(s) [²⁷]. Various weights can return spurious optimum points, which could be far from genuine. In addition, if the objective function is in the form of a multi-objective function, then it may return a Pareto optimal front issue, which is difficult to resolve and refine in higher-dimensional space [¹³]. The optimization codes must minimize the objective function by updating the input parameters either randomly or by using gradient-based schemes. The sensitivity of the eigenvalues and eigenvectors (features) of any system to unknown parameters can be expressed as follows [¹²]:

(4)\frac{\partial \langle {\left(2\pi f_i\right)}^2\rangle }{\partial p_j}={\left\{{\varphi }_i\right\}}^{\mathrm{T}}\left(\frac{\partial \left[K\right]}{\partial p_j}-\langle {\left(2\pi f_i\right)}^2\rangle \frac{\partial \left[M\right]}{\partial p_j}\right)\left\{{\varphi }_i\right\},

(5)\frac{\partial \left\{{\varphi }_i\right\}}{\partial p_j}=-{\displaystyle \sum _{s=1}^N\left(\frac{{\left\{{\varphi }_s\right\}}^{\mathrm{T}}\left(\frac{\partial \left[K\right]}{\partial p_j}-\langle {\left(2\pi f_i\right)}^2\rangle \frac{\partial \left[M\right]}{\partial p_j}\right)\left\{{\varphi }_i\right\}}{{\left(2\pi f_s\right)}^2-{\left(2\pi f_i\right)}^2}\right)\left\{{\varphi }_{\mathrm{s}}\right\}},

where [M], [K], p_j, f_i, and j_i are the mass matrix, stiffness matrix, jth unknown parameter, ith NF, and mode shape of the system, respectively. As can be seen from Eqs. (4) and (5), any change in the stiffness and mass matrices of the system can lead to changes in the system features. In addition, as the relationships between the system parameters and their features are nonlinear, the degree of nonlinearity can be estimated using numerical simulations and regression analyses.

Proposed methodology

The proposed methodology is designed to identify the minimum number of essential global and local features to determine the locations of the measurement node and uncertain system parameters. A flowchart of the proposed methodology is presented in Fig. 1. The methodology is divided into a multi-stage process as follows.

STAGE (1): construct an EVP with a PDE/FEM

The first step is to construct an EVP for the desired structural system, which can be solved analytically or numerically. For simple discrete and continuous elastic media such as rods, beams, plates, or shells with simple BCs, it is possible to analytically solve the EVP. In these cases, the behavior of the system can be expressed by a partial differential equation (PDE), and then a nonlinear characteristic equation must be derived for further input-output (I/O) analysis. Alternatively, a detailed FEM of the structure can be constructed to numerically solve its EVP.

STAGE (2): design of experiments

As the relationship between the system’s features and parameters is implicit, design of experiments (DOE) can be exploited to generate several input parameters. These trial points are various combinations of the system’s uncertain parameters within their FR [²⁸]. The number of samples depends on the number of design variables or RVs, the size of the FR, the number of DOFs of the boundary variables, and the acceptable cross correlation between various RVs. In this research, Latin hypercube sampling (LHS) was used for the DOE process. This method shuffles samples of multi-dimensional variables randomly within the FR without any priority [²⁹]. With this setup, each variable is distributed uniformly, and the least correlation exists between any pair of variables to cover the entire FR appropriately.

STAGE (3): solve and save EVPs

In this step, the constructed EVP must be solved for all DOEs. The physical parameters of the entire system, \underset{\_}{X}=\{x_{\mathrm{1}}\dots x_M\}, construct a vector of M number of unknowns. Each output of the system y_n is usually either an eigenvalue or an eigenfunction. The outputs of the entire system, \underline{Y}=\{f_i,{\varphi }_i:i=\mathrm{1}:N\}, construct 2N system features. It has been strongly suggested that, to determine a viable solution, the number of selected features must be equal to or greater than the number of parameters that need to be identified [²⁰]. In other words, in order to estimate M system parameters authentically, at least N modes must be excited, where (N≥M). The input parameters and corresponding output global features and local features are used to construct regression surfaces hereafter.

STAGE (4): find sensitive nodes

In this step, collection points for local information are identified. These are points where spots/nodes oscillate the most when the parameters of the system change. The local information depends on the location as well as the directions (x, y, and z) of motion at the measurement points. Because it is computationally expensive to collect all the mode shapes of a system, it would be more practical to simply find the most sensitive nodal information [³⁰]. The mode shape component, j_ij, represents the amplitude of the ith mode at the jth DOF. The selection of the node number and its DOF in each mode depends on its sensitivity to various system parameter variabilities. In this section, we present the modal transmissibility (MT) to represent the local feature hereafter. The MT is defined by the ratio of the mode shape components of the jth and kth DOFs in the ith mode as follows:

(6)M{T}_i^{jk}={\displaystyle \frac{{\varphi }_{ij}}{{\mathrm{phi;}}_{ik}}}\mathrm{.}

In practice, this local variable can be estimated by power spectral density transmissibility (PSDT). First, the ratio of the PSD between the time domain signals of y_j(t) and y_k(t) with a reference signal y_r(t) is defined by the following formula [³¹,³²]:

(7)T_{jk}^r(f)={\displaystyle \frac{G_{jr}(f)}{G_{kr}(f)}},

where G_jr(f) and G_kr(f) are the cross PSDs of the corresponding signals. The PSDT is a frequency-dependent function. Around the ith NF of the system, this function returns the MT_i as follows [³³–³⁵]:

(8)\underset{f\to f_i}{\lim }\hspace{0.17em}{T}_{jk}^r(f=f_i)={\displaystyle \frac{{\varphi }_{ij}}{{\varphi }_{ik}}}=M{T}_i^{jk}\mathrm{.}

Numerically, the entire mode shape can be easily collected for further statistical processing. Then, after finding the best sensitive nodes on the structure, the data at these nodes can be mounted for the feature extraction process during sensor collection. The minima and maxima of mode shapes at any node are associated with the smallest and largest energies, respectively. If the minimum of a mode shape at a node is close to zero, that node attains the least energy and is called a fixed/modal node. On the other hand, if the maximum of the mode shapes is an antinode, this means that the highest energy is concentrated on them during excitation [²⁶,³⁶]. It was found that the most appropriate locations for data collection were the points that attained the highest standard deviation (STD) while being far away from the fixed points or minimum. A high STD means that the nodes oscillate the most when the system parameters change and attain the highest sensitivity to the predefined system parameters. The source of a high STD of any node could be at least one uncertain parameter because the system parameters vary simultaneously.

STAGE (5): create global and local information RSs

Because the FEM is a type of black box and the interrelationship of various I/O parameters is unknown, the response surface (RS) method is presented in this work as a substitute for the FEM [²⁸,³⁷]. It will be shown that highly nonlinear relationships exist between the local and global features of the system and system parameters. Owing to the derivability and tractability provided by polynomial expressions, polynomials are used in the curve-fitting process. It has been proven that the number of regression coefficients required for precise curve fitting grows exponentially as the polynomial degree increases [⁹,¹¹]. The general high-order polynomial RS model plus its interaction terms for the system feature, y can be expressed as [³⁸,³⁹]

(9)y\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_M\right)={\displaystyle \sum }_{i_{\mathrm{1}}=\mathrm{0}}^p{\displaystyle \sum }_{i_{\mathrm{2}}=\mathrm{0}}^{p-i_{\mathrm{1}}}{\displaystyle \sum }_{i_{\mathrm{3}}=\mathrm{0}}^{p-i_{\mathrm{1}}-i_{\mathrm{2}}}\cdots {\displaystyle \sum }_{i_M=\mathrm{0}}^{p-i_{\mathrm{1}}-i_{\mathrm{2}}\cdots i_{M-\mathrm{1}}}

where{\beta }_{i_{\mathrm{1}}{i}_{\mathrm{2}}\dots i_M}x_n, e, M, and p are the regression coefficients, nth independent variable or predictor parameter, Gaussian error term, total number of independent variables, and highest order of polynomials, respectively. There are two major indicators of the goodness and accuracy of curve-fitted surfaces/functions. The coefficient of determination, R², is the first indicator and represents the ratio of the sum of the square regression (SSR) to the total sum of the square (SST) as follows [¹]:

(10)R^{\mathrm{2}}={\displaystyle \frac{SSR}{SST}}={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n{\left({\stackrel{�}{y}}_i-\overline{y}\right)}^{\mathrm{2}}}}{{\displaystyle \sum _{i=\mathrm{1}}^n{\left(y_i-\overline{y}\right)}^{\mathrm{2}}}}}\mathrm{.}

This coefficient (R²) can vary from 0.0, which indicates no accuracy, to 1.0. The second indicator is the root mean square error (RMSE), which evaluates the overall accuracy of the fitting surface and is defined as follows:

(11)RMSE=\sqrt{{\displaystyle \frac{SST}{n}}}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{n}}}\mathrm{.}

Usually, the RS method is executed by a second-order (p = 2, quadratic) or at most third-order (p = 3, cubic) polynomial, owing to the simplicity in estimating its coefficients [⁴,⁴⁰]. However, when considering several uncertain parameters simultaneously, satisfactory R² can barely be achieved by a low polynomial order. However, the real degree of nonlinearity between the system features and its uncertain parameters is demanded in this proposed method. Practically, it is impossible to achieve an exact RS to express the I/O relationship. Therefore, the order of the polynomials was gradually increased until a high accuracy (R²≥0.99) and minimum RMSE were achieved. The overall problem definition can be expressed as follows: If there are M system parameters (x_i) that need to be identified, N features, (y_j) must be extracted to correctly identify (x_i). In other words, the system parameter vector is \underset{\_}{X}=\{x_{\mathrm{1}}\dots x_M\}, and the system feature vector is \underset{\_}{Y}=\{y_{\mathrm{1}}\dots y_N\}. The direct and inverse problems, respectively, are as follows:

(12)\underline{X}=\{x_{\mathrm{1}}\dots x_M\}\stackrel{\mathrm{y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\mathrm{s}}{\to }\underline{Y}=\{y_{\mathrm{1}}\dots y_N\},

(13)\underline{Y}=\{y_{\mathrm{1}}\dots y_N\}\stackrel{\mathrm{y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\mathrm{s}}{\to }\underline{X}=\{x_{\mathrm{1}}\dots x_M\}\mathrm{.}

In the direct approach, either the PDE or FEM returns the system features based on predefined model parameters, while the inverse problems aim to find system parameters from extracted model features. The crucial issue in almost all inverse problem solvers is determining whether the returned result is a genuine or a false one, especially for a high-dimensional domain in which there is no way to investigate the optimum domain.

STAGE (6): create a myriad discrete set

To check the uniqueness of the FR, the infinite continuous multidimensional space is discretized to a plentiful number (N_D) of random trial points that cover the FR accurately using the same method as in Stage 2 while utilizing the highest number of trial points that the computer can handle for searching. Because the parallel coordinate can only plot a discrete set of multidimensional points, this step must be executed for further visualization. As can be seen in the second and third examples, N_D = 1E6 random trial points were generated by the DOE process for the discretization of the domains.

STAGE (7): I/O evaluation of a myriad set from RSs

In this stage, the myriad set must be evaluated by the RSs estimated from the global and local information (Stage 5). This process is computationally trivial because there is no longer a need to solve EVPs; instead, one just needs to evaluate RSs. At the end of this step, the I/O values of the myriad points can be estimated in a metadata format\{\underline{\mathit{X}}:\underline{\mathit{Y}}\}. Because the continuous FR was discretized in a multidimensional space, there is still a gap between each sequential input \underset{\_}{X} and output \underset{\_}{Y} of the system. In this proposed methodology, instead of putting ambiguous weights, such as those in Eq. (1), into the system features, each feature is assigned a unique weight based on its own range in the FR. Because each feature varies independently of other features as the input parameter changes in the FR, selection of the desired points is based on an equal percentage of range/variability for each feature. In the constructed set \{\underline{\mathit{X}}:\underline{\mathit{Y}}\}, the feature gap percentage (FGP) is presented as the quotient of the maximum difference between two successive or sequentially ordered features in the set divided by its own potential range. This can be considered the existent discretization gap due to the sampling at the previous step. The scalar FGP index for the constructed set can be evaluated by

(14)FGP=Max\{Distance[Ordered(y_i)]/Range(y_i);i=\mathrm{1}:N\}\mathrm{.}

The FGP can be considered a percentage by which each feature in the constructed set can increase or decrease to its neighborhood point. It must be mentioned that when the number of NDs increases, the FGP decreases by the same amount. The maximum of the individual feature gap must be chosen to ensure that the search process will find at least one sample.

STAGE (8): find minimum required modal information

This stage is the most critical; it indicates the number of essential eigen solutions required to identify the system parameters. It is always desirable to identify system parameters with the smallest number of excited modes because it is costly to excite many modes of a large structure. The additional extracted features could be redundant and could ruin the results because there is less accuracy in the higher mode features in practice [³⁰]. If we know how many features are sufficient, then more energy must be used in extracting them more precisely based on the required accuracy rather than finding redundant ones. In Stage 5, it was found that the relationships between system parameters and their features are highly nonlinear. It, therefore, needs to be determined whether a unique intersection domain of local and global features can be achieved. For two-dimensional (2D) functions, a contour plot can visually show the intersection of the local and global features, as in the first example. Whereas if a feature be function of more than two independent variables, it is challenging to plot that function in common 3D Cartesian domain. A novel contribution of this work is the use of the parallel coordinate for SMV purposes, and its feasibility is demonstrated in the second and third examples. The parallel coordinate is a well-known geometric toolbox for the visualization of discretized high-dimensional data [⁴¹–⁴³].

The eighth stage comprises several steps, as shown in the dotted box at the lower end of the flowchart in Fig. 1. The goal of this stage is to find the minimum number of modes Nm that can provide sufficient constraints to authentically find M system parameters. There are several reasons why every conventional optimization code returns several optimum points. First, based on the definition of the conventional objective function described by Eq. (1), different weights return different parameters. Second, as extracted features are obtained repeatedly during experiments, it is inevitable that these features will have some variations even in the same environment. Finally, there is no way to determine whether these returned optimum points are genuine. Therefore, it is anticipated that the proposed parallel coordinate will find a unique domain (UD), which means that each returned system parameter is stable around a single neighborhood. As mentioned earlier, all optimization codes provide several optimum points that could vary in FR, but the UD returns a unique intersection around single system parameter values.

Consequently, an iterative loop was defined to reach a UD from the smallest number of modes Nm= 1, as shown in Fig. 1. In each iteration, a trial point of system parameters \underline{\mathit{X}}(\mathit{i}) was randomly generated within the FR in the same manner as in the second stage. Then, the system features \underline{\mathit{Y}}(\mathit{i}) that contain only Nm modal information could be estimated in a similar way as in the third stage. Now, even if the exact value of \underline{\mathit{Y}}(\mathit{i}) does not exist in the constructed discrete myriad set \{\underline{\mathit{X}}:\underline{\mathit{Y}}\}, the algorithm will still be able to find the optimum system parameters \underline{{\mathit{X}}_o} in the set as follows:

(15)OptimumPoints=\mathrm{Find}\left\{\underline{{\mathit{X}}_o}|\underline{\mathit{Y}}(\mathit{i})-\underline{B_Y}\le \underline{\mathit{Y}}\le \underline{\mathit{Y}}(i)+\underline{B_Y}\right\},

where the feature bound vector (\underset{\_}{B_Y}) is defined to provide a fine interval:

(16)\underline{B_Y}=FGP*Range\left(\underline{\mathit{Y}}\right)\mathrm{.}

The FGP gives weights to the system features by assigning a degree of importance to each feature based on its variability inside FR. Afterwards, Eqs. (15)–16 return several optimum points \underline{{\mathit{X}}_o} from the myriad set. Then, the algorithm increases the number of modes Nm during the search process of the defined loop automatically until a UD of \underline{{\mathit{X}}_o} is achieved. The UD can be defined when the distance between various \underline{{\mathit{X}}_o} becomes less than an acceptable value. The UD for the first example is displayed with a counter plot as two independent variables exist. For the second and third examples, the UD is investigated by parallel coordinates because more than two independent variables exist in the RSs.

Examples and verifications

This section aims to thoroughly demonstrate all the required steps based on the proposed methodology depicted in Fig. 1. Three examples, ranging from a simple structural model to more complex ones, were chosen. The first and second examples present prismatic beams supported on various BCs described by the Euler-Bernoulli theory, and the third example is a full three-dimensional FEM of a highway bridge superstructure. In this research, MATLAB® was adopted for the DOE process, root-finding of nonlinear characteristic equations, high-order polynomial curve-fitting, and data visualization. In the third example, ABAQUS and the Abaqus2Matlab toolbox [⁴⁴] were exploited for the eigen solver machine and data collection purposes.

Example 1: prismatic Euler-Bernoulli beam plus rotational spring supports at each end (two system parameters)

A well-known EVP is a continuous prismatic beam of a homogeneous material described by modulus of elasticity E, second moment of area I, density r, length L, and cross section A. The flexural displacement of the Euler-Bernoulli beam can be written as the following PDE [⁴⁵,⁴⁶]:

(17)\rho A{\displaystyle \frac{{\partial }^{\mathrm{2}}w}{\partial t^{\mathrm{2}}}}+EI{\displaystyle \frac{{\partial }^{\mathrm{4}}w}{\partial x^{\mathrm{4}}}}=\mathrm{0.}

Using the separation of variables technique, this PDE has the following solution:

(18)w(x,t)=\varphi (x).q(t)={\displaystyle \sum }_{n=\mathrm{1}}^{\infty }{\varphi }_n(x).q_n(t)={\displaystyle \sum }_{n=\mathrm{1}}^{\infty }{\varphi }_n(x).e^{i(\mathrm{2}\pi f_n)t},

where j_n, q_n, and f_n are the nth mode shape, frequency equation, and resonant frequency, respectively. Substituting Eq. (18) into Eq. (17) and making a few calculation simplifications yields the following univariable equations:

(19){\displaystyle \frac{d^{\mathrm{2}}{q}_n(t)}{\partial t^{\mathrm{2}}}}+q_n(t)=\mathrm{0},

(20){\displaystyle \frac{d^{\mathrm{4}}{\varphi }_n(x)}{\partial x^{\mathrm{4}}}}-{\beta }_n^{\mathrm{4}}{\varphi }_n(x)=\mathrm{0},

where β_n is the nth modal coefficient that depends on the nth frequency coefficient (modal number) λ_n. This modal coefficient is related to the material properties as well as the geometry of the beam as follows:

(21){\beta }_n^{\mathrm{4}}={\left({\displaystyle \frac{{\lambda }_n}{L}}\right)}^{\mathrm{4}}={(\mathrm{2}\pi f_n)}^{\mathrm{2}}{\displaystyle \frac{\rho A}{EI}},

(22)f_n={\displaystyle \frac{{\lambda }_n^{\mathrm{2}}}{\mathrm{2}\pi L^{\mathrm{2}}}}\sqrt{{\displaystyle \frac{EI}{\rho A}}}\mathrm{.}

The general solution of the EVP described by Eq. (20), that is, the mode shapes of the system, can be calculated as follows [⁴⁷]:

(23){\varphi }_n(x)=A_n.\sin ({\beta }_{n}x)+B_n.\cos ({\beta }_{n}x)+C_n.\sinh ({\beta }_{n}x)+D_n.\cosh ({\beta }_{n}x),

where A_n, B_n, C_n, and D_n are the constant coefficients for each mode, which can be estimated based on various BCs. By substituting various BCs into the latter equation, a nonlinear univariable equation, namely the characteristic equation with infinite roots, is obtained. Its roots are modal numbers λ_n, which are related to the system NF by Eq. (21).

The first example, which is a Euler-Bernoulli beam with two rotational springs, is shown in Fig. 2(a). The four BCs of this model are as follows [⁴⁸]:

(24)\{\begin{array}{c}{\phi }_n\left(x=0\right)=0\hfill \\ \hfill \\ {\phi }_n\left(x=L\right)=0\hfill \\ \hfill \\ {\displaystyle \frac{{\mathrm{d}}^2{\phi }_n\left(x=0\right)}{\mathrm{d}{x}^2}}={\displaystyle \frac{K_{\theta \mathrm{l}}}{EI}}{\displaystyle \frac{\mathrm{d}{\phi }_n\left(x=0\right)}{\mathrm{d}x}}\hfill \\ \hfill \\ {\displaystyle \frac{{\mathrm{d}}^2{\phi }_n\left(x=L\right)}{\mathrm{d}{x}^2}}={\displaystyle \frac{-K_{\theta \mathrm{r}}}{EI}}{\displaystyle \frac{\mathrm{d}{\phi }_n\left(x=L\right)}{\mathrm{d}x}}\hfill \end{array}\mathrm{.}

Substituting the BCs shown in Eq. (24) into the general solution of Eq. (23) and making several parametric manipulations yields the characteristic equations (frequency equation):

(25)\begin{array}{l}{K}_{\theta l}^{*\mathrm{2}}+K_{\theta l}^{*}{\displaystyle \frac{{\lambda }_n[\mathrm{1}+\alpha ][\sin ({\lambda }_n)\cosh ({\lambda }_n)-\cos ({\lambda }_n)\sinh ({\lambda }_n)]}{\alpha [\mathrm{1}-\cos ({\lambda }_n)\cosh ({\lambda }_n)]}}+\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\mathrm{2}{\lambda }_n^{\mathrm{2}}\sin ({\lambda }_n)\sinh ({\lambda }_n)}{\alpha [\mathrm{1}-\cos ({\lambda }_n)\cosh ({\lambda }_n)]}}=\mathrm{0}\end{array}\mathrm{.}

This equation is remarkably nonlinear, and its roots need to be evaluated for various spring coefficients using a robust root-finding toolbox. Then, by substituting the roots of this equation as well as the modal unknown coefficients into Eq. (23), the normalized mode shape can be determined as follows [⁴⁹]:

(26)\begin{array}{l}{\mathrm{phi;}}_n(\xi )=\sin ({\lambda }_n\xi )-\sinh ({\lambda }_n\xi )+\left({\displaystyle \frac{\sinh ({\lambda }_n)-\sin ({\lambda }_n)}{\cos ({\lambda }_n)-\cosh ({\lambda }_n)-{\displaystyle \frac{\mathrm{2}{\lambda }_n}{K_{\theta l}^{*}}}\sinh ({\lambda }_n)}}\right)\times \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left[\cos ({\lambda }_n\xi )-\cosh ({\lambda }_n\xi )-{\displaystyle \frac{\mathrm{2}{\lambda }_n}{K_{\theta l}^{*}}}\sinh ({\lambda }_n\xi )\right]\end{array},

where \alpha ={\displaystyle \frac{K_{\theta r}}{K_{\theta l}}}, K_{\theta l}^{*}={\displaystyle \frac{K_{\theta l}}{(EI/L)}}, K_{\theta r}^{*}={\displaystyle \frac{K_{\theta r}}{(EI/L)}}, and \xi =x/L are the rotational spring ratio, dimensionless left-side rotational spring, dimensionless right-side rotational spring, and dimensionless beam length, respectively. As Eq. (25) and Eq. (26) are evaluated by solving a pre-defined PDE and then become dimensionless, they are applicable to all similar prismatic beams with the same BCs. In this example, two dimensionless spring constants (K_{\theta \mathrm{l}}^{*}, K_{\theta \mathrm{r}}^{*}) are considered as the system parameters. Any combination of these dimensionless constants in the characteristic Eq. (25) will yield infinite modal parameters λ_n as the equation roots. Then, those roots will construct infinite mode shapes of the system by Eq. (26). As only NFs could be extracted experimentally and the characteristic equation roots are the modal numbers, the following formula, which relates modal numbers λ_n to the normalized NFs, was used:

(27){\displaystyle \frac{f_n}{f_{\mathrm{1}}}}={\displaystyle \frac{{\lambda }_n^{\mathrm{2}}}{{\lambda }_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{.}

This formula can be confirmed by Eq. (22); somehow, all NFs are normalized to the first NF. Therefore, numerical as well as experimental results could be completely related. In Fig. 3(a), the correlation plot of the 500 samples generated by the LHS algorithm is shown.

By substituting these samples into Eq. (25), a collection of modal numbers λ_n based on various (K_{\theta \mathrm{l}}^{*}, K_{\theta \mathrm{r}}^{*}) was estimated numerically by the highest precisions the root finding toolbox could handle. Then, the normalized mode shapes were evaluated based on Eq. (26); the overlaid mode shapes of the first three modes and their statistical information are plotted in Fig. 4. For this example, nodes at ζ= 0.2 and ζ= 0.8 were chosen for the local information collection (MT) based on Eq. (6), as these nodes attain relatively high STD while being located far away from the nodal points of the first three modes. The local information in this example, i.e., MT1, MT2, and MT3, are the MT features at those nodes in the first, second, and third modes, respectively. For RS fitting, global and local information must be chosen; they intersect at one region/domain. The RS of a few initial global and local features of the samples as well as curve-fitted data are shown in Fig. 5. In this example, (p= 4) was applied as the highest polynomial order to obtain a high coefficient of determination (R2≥0.99). As there are only two independent parameters, the intersection domain can be visually identified.

The contour plots of the feature RSs by specific trial values are shown in Fig. 6(a) and Fig. 6(b) for FGP = 0.01 and FGP = 0.05, respectively. From the Fig. 6 contour plots, three important conclusions can be made. First, both local and global information are required to update and identify system parameters at a UD. In other words, neither global nor local information individually can specify a UD; however, the region where the two intersect can. Second, as the feature value uncertainty FGP increases, a broader domain of optimal points can be returned by the UD. This plot shows that, by enlarging the feature uncertainty, the provided UD swells and consequently proliferates the potential optimum points. Finally, it is clear from the contour plot that at least two modes must be excited and that their global and sensitive local features are able to return the correct system parameters (N_m = 2). The other modal information is redundant, as shown in Fig. 6.

Example 2: prismatic Euler-Bernoulli beam plus translational and rotational spring supports at each end (four system parameters)

In this example, both translational and rotational springs at the bearing sides are considered simultaneously as the most general case for the Euler-Bernoulli beam BCs, as shown in Fig. 2(b). The four BCs in this problem are as follows [⁴⁸]:

(28)\{\begin{array}{l}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{3}}{\phi }_n\left(x=\mathrm{0}\right)}{\mathrm{d}{x}^{\mathrm{3}}}}=-{\displaystyle \frac{K_{X\mathrm{l}}}{EI}}{\phi }_n\left(x=\mathrm{0}\right)\hfill \\ \hfill \\ {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{\phi }_n\left(x=\mathrm{0}\right)}{\mathrm{d}{x}^{\mathrm{2}}}}={\displaystyle \frac{K_{\theta \mathrm{l}}}{EI}}{\displaystyle \frac{\mathrm{d}{\phi }_n\left(x=\mathrm{0}\right)}{\mathrm{d}x}}\hfill \\ \hfill \\ {\displaystyle \frac{{\mathrm{d}}^{\mathrm{3}}{\phi }_n\left(x=L\right)}{\mathrm{d}{x}^{\mathrm{3}}}}={\displaystyle \frac{K_{X\mathrm{r}}}{EI}}{\phi }_n\left(x=L\right)\hfill \\ \hfill \\ {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{\phi }_n\left(x=L\right)}{\mathrm{d}{x}^{\mathrm{2}}}}={\displaystyle \frac{-K_{\theta \mathrm{r}}}{EI}}{\displaystyle \frac{\mathrm{d}{\phi }_n\left(x=L\right)}{\mathrm{d}x}}\hfill \end{array}\mathrm{.}

Substituting the BCs shown in Eq. (28) into the general solution of Eq. (23) and making a few parametric manipulations gives the following characteristic Eqs [⁴⁹]:

(29)\begin{array}{l}{\lambda }_n^{\mathrm{8}}[\mathrm{1}-\cos ({\lambda }_n)\cosh ({\lambda }_n)]-{\lambda }_n^{\mathrm{7}}[(K_{\theta l}^{*}+K_{\theta r}^{*})(\sin ({\lambda }_n)\cosh ({\lambda }_n)+\cos ({\lambda }_n)\sinh ({\lambda }_n))]\\ -{\lambda }_n^{\mathrm{6}}[\mathrm{2}{K}_{\theta l}^{*}{K}_{\theta r}^{*}\sin ({\lambda }_n)\sinh ({\lambda }_n)]+{\lambda }_n^{\mathrm{5}}[(K_{\mathrm{Xl}}^{*}+K_{\mathrm{Xr}}^{*})(\cos ({\lambda }_n)\sinh ({\lambda }_n)-\sin ({\lambda }_n)\cosh ({\lambda }_n))]\\ +{\lambda }_n^{\mathrm{4}}[(K_{\theta l}^{*}{K}_{\mathrm{Xl}}^{*}+K_{\mathrm{Xr}}^{*}{K}_{\theta r}^{*})(\mathrm{1}+\cos ({\lambda }_n)\cosh ({\lambda }_n))+\mathrm{2}(K_{\theta r}^{*}{K}_{\mathrm{Xl}}^{*}+K_{\mathrm{Xr}}^{*}{K}_{\theta l}^{*})(\cos ({\lambda }_n)\cosh ({\lambda }_n))]\\ +{\lambda }_n^{\mathrm{3}}[K_{\theta l}^{*}{K}_{\theta r}^{*}(K_{\mathrm{Xl}}^{*}+K_{\mathrm{Xr}}^{*})(\cos ({\lambda }_n)\sinh ({\lambda }_n)+\sin ({\lambda }_n)\cosh ({\lambda }_n))]+\\ \mathrm{2}{\lambda }_n^{\mathrm{2}}[K_{\mathrm{Xl}}^{*}{K}_{\mathrm{Xr}}^{*}\sinh ({\lambda }_n)\sin ({\lambda }_n)]+{\lambda }_n[K_{\mathrm{Xl}}^{*}{K}_{\mathrm{Xr}}^{*}(K_{\theta l}^{*}+K_{\theta r}^{*})(\sin ({\lambda }_n)\cosh ({\lambda }_n)-\cos ({\lambda }_n)\sinh ({\lambda }_n))]\\ +K_{\theta l}^{*}{K}_{\theta r}^{*}{K}_{\mathrm{Xl}}^{*}{K}_{\mathrm{Xr}}^{*}[\mathrm{1}-\cos ({\lambda }_n)\cosh ({\lambda }_n)]=\mathrm{0}\end{array}

This equation is nonlinear, and its roots depend on K_{\mathrm{Xl}}^{*}={\displaystyle \frac{K_{\mathrm{Xl}}}{EI/L^{\mathrm{3}}}} and K_{\mathrm{Xr}}^{*}={\displaystyle \frac{K_{\mathrm{Xr}}}{EI/L^{\mathrm{3}}}}, which are the left and right translational normalized spring coefficients, respectively. The two rotational spring coefficients K_{\theta l}^{*} and K_{\theta r}^{*}are the same as in example 1. In this example, four dimensionless spring constants (K_{\mathrm{Xl}}^{*}, K_{\mathrm{Xr}}^{*}, K_{\theta l}^{*},K_{\theta r}^{*}) are considered as the system parameters. The mode shapes are determined by inserting the roots of Eq. (29) into the mode-shape formula. Because it is a very long equation, it will not be shown in the paper. However, it is available in the literature [⁵⁰,⁵¹].

In Fig. 3(b), the correlation plot of the 1000 samples generated by the LHS algorithm is shown. The mode shapes of the system based on various sample values can be evaluated. The overlaid mode shapes of the first five modes and their statistical information are plotted in Fig. 7. In this example, nodes at ζ= 0.1 and ζ= 0.9 were chosen for local information collection (MT) based on Eq. (6) because they met the high STD criteria while being far away from the nodal points of the initial five modes of the beam. The local information in this example, MT1, MT2, MT3, MT4, and MT5, are the MT features at those nodes in the first through fifth modes, respectively. In this example, (p = 5) was applied as the highest polynomial order to obtain a high coefficient of determination (R²≥0.99). The order of the polynomials is one degree more than that in the previous example, where rotational springs were considered. Although the RSs cannot be plotted as four independent variables, the statistical information for the global and local features of the 1000 samples is shown in Fig. 8 to check their sensitivity. The intersection domain must be investigated using the parallel coordinate toolbox of the myriad data set to find the UD. To find the essential information for correct parameter identification, the myriad discretization step was performed with ND= 1E6 points provided by LHS and FGP= 0.01. In Fig. 9, parallel coordinate plots are shown for the myriad and the gradually increasing Nm from just two modes up to four modes of trial features. Several different trial feature points were used to ensure that there was enough information to find the UD. Based on this toolbox, the four modes must be excited and their global and local information at the pre-mentioned nodes must be collected in order to find the UD (Nm= 4). In this case, the information in the fifth mode is redundant and can be disregarded during the experimental feature extraction process.

Example 3: highway bridge super structure plus rotational BCs (six system parameters)

This structural model was selected as a benchmark for a more complex structure. It is one span of a highway bridge constructed by a concrete deck on four girders and 12 cross beams, as shown in Fig. 2(c). The generalized EVP can be expressed by the following Eq [⁵²]:

(30)[[K]-{(\mathrm{2}\pi f_i)}^{\mathrm{2}}[M]]\{{\varphi }_i\}=\mathrm{0},

where M, K, F_i, and j_i are the mass matrix, stiffness matrix, ith NF, and mode shape of the system, respectively.

In this example, six uncertain parameters were considered for the SMV process, \{E_{\mathrm{s}},E_{\mathrm{c}},{\rho }_{\mathrm{s}},{\rho }_{\mathrm{c}},K_{\theta \mathrm{l}},K_{\theta \mathrm{r}}\}, representing the Young’s modulus of elasticity of stainless steel used in the girders and cross-beams, Young’s modulus of elasticity of the concrete used in the deck, density of stainless steel, density of concrete, left side of the rotational stiffness of the bearing system, and right side of the rotational stiffness of the bearing system, respectively. The three-dimensional FEM of the bridge was simulated by 7929 elements of the 4-node quadrilateral elastic shell element in Abaqus®. The overall bridge model had 9028 nodes after it was meshed in the software. Because the FEM is deployed as the EVP solver, the number of elements was obtained by starting from an initial guess and proceeding to a finer mesh until the relative variation of the required NFs became less than the acceptable tolerance defined by the analyst.

In this study, the mesh sensitivity was determined from the initial guess of the mesh size until the six initial NF variations were less than 0.1% of the previous iteration. In Fig. 3(c), the correlation plot of the 500 samples generated by the LHS algorithm is shown. Then, all DOE combinations were inserted as inputs to the FEM in Abaqus®, and all the eigen solutions were collected automatically for further analysis. The six initial bending mode shapes and their STD contour plots are shown in Fig. 10. Only bending modes were analyzed because these modes are more sensitive to the system parameters. The DOFs required for estimating MT_i were chosen at ζ= 0.135 and ζ= 0.865 at the middle longitudinal line of the deck. In Fig. 11, the statistical information of the selected global and local features is shown and is used to investigate the sensitivity and precision of each feature. In this example, (p = 4) was applied as the highest polynomial order to obtain a high coefficient of determination (R²≥0.99). The intersection domain must be investigated using the parallel coordinate toolbox of the myriad set to find the UD. To determine the essential information for correct parameter identification, the myriad discretization step was performed with N_D = 1E6 provided by LHS and FGP = 0.01. In Fig. 12, the parallel coordinate plot of the myriad data set as well as the gradually increasing constraints (Nm) from four to six modes of a trial feature are plotted. Again, a few various trial feature points were used to determine whether the number would provide sufficient information to find the UD. Through the parallel coordinate toolbox, it was found that at least six modes must be excited, and then their global and local information at the pre-mentioned nodes must be collected to authentically find the UD (Nm = 6).

In this example, the standardization process was also performed for better visualization using parallel coordinates. Because the system parameters have different units and ranges, each system parameter in the myriad set must be standardized by (Z_i={\displaystyle \frac{X_i-{\mu }_{X_i}}{{\sigma }_{X_i}}}). There are a few notes that an analyst must keep in mind when using any FEM software during the required EVP process, such as element type, mesh size, and BCs, in order to minimize discretization errors. If inappropriate values are used for any of these simulation sources can mislead the optimization algorithm toward spurious points.

Possible experimental implementation

This section presents the experimental application of the theoretical/numerical algorithm presented in Fig. 1. As indicated in Section 2, the algorithm presented in Fig. 1 provides the locations of the critical nodes that have high energies and that are far away from the nodal points. While the algorithm can identify more candidate nodes if needed, the three examples solved in Section 3 showed that even two nodes were sufficient to identify six system parameters. The algorithm also creates an RS between each modal feature and the system parameters for large cases within the FR. Having this information available will assist the analyst in conducting experiments and identifying the unknown system parameters as accurately as possible. During the experiments, the acceleration data from the sensors were used to identify the modal features of the system represented by NFs as the global information and MTs as the sensitive local information [²⁶,⁵³]. The NFs can be identified from the PSD peaks of the signals, and the ratio between the PSD magnitudes of the two sensors at the different peaks (NFs) will be used to calculate the MT at each NF [³⁵,⁵⁴].

Although the mode shapes for all DOEs must be evaluated numerically to be used in the statistical analysis, only specific mode shape indices are experimentally extracted to construct in situ MTs. The computational cost of mode shape extraction by computer-based FE code is negligible compared to the experimental cost associated with mounting many sensors. Rather, a few sensors can be deployed to extract similar structural information (NFs and MTs). Once the real system features are determined, the next step is to determine the system parameters using an inverse problem. However, this process does not require solving a computationally expensive inverse problem; instead, it will rely on matching the set of experimental features with the myriad data set in the RSs. An iterative process is performed until the experimental features match the closest system parameters presented in the myriad data set. Because of the discrete nature of the myriad data set, there is no guarantee that an exact matching solution between the myriad random set and the experimental features exist [⁵⁵]. Nevertheless, the bound vector (\underset{\_}{B_Y}), which is based on the statistical range of the features in FR as defined by Stage 8 of the algorithm, will be used to exquisitely perturb the magnitude of the experimental features until it matches the closet parameters in the myriad set. Eventually, the returned results can be visually inspected by the parallel coordinate tool [⁵⁶].

Discussion and conclusions

This paper presents a novel approach for performing a structural model updating process that minimizes system uncertainties and determines a unique solution for the system’s parameters and BCs using a minimum number of measurement points and system features. A new approach for assigning objective weights to each term in a multi-objective function is proposed based on the system’s local and global information, and its sensitivity and range in the FR. Assigning objective weights to each term will swiftly increase the efficiency of the optimization process. Second, in contrast to conventional SMV methods, which select arbitrary positions for data collection, the proposed method used a new approach to identify the most sensitive nodes for sensor placement. Third, the range-based weights and the parallel coordinate toolbox were exploited as two unprecedented approaches for SMV applications to identify an authentic UD with a minimum number of required features. Finally, in comparison with existing methods, high-order polynomial RSs are suggested to measure valid degrees of nonlinearity between system parameters and their features. These nonlinearities are the reason that most updating methodologies are unable to find genuine optimum points.

In the first and third examples, where just rotational BCs were considered as the system parameters, fourth-order polynomials were found to be appropriate for estimating acceptable RSs. However, in the second example, fifth-order polynomials were required to obtain a similar accuracy on the optimum point. The interaction between both rotational and translational DOFs is the main reason why higher-order polynomials were utilized to find UD of optimum points with acceptable accuracy. These RSs with high accuracy are deployed as pattern recognition functions to evaluate myriad random points, instead of running finite element (FE) code for enormous lengths of time. However, a caveat exists when using high-order polynomials to fit any data set. Extreme high-order polynomials can cause severe ringing between the data points and further spurious results. Therefore, the algorithm must gradually enlarge the order of the polynomial from low to high until an acceptable criterion is achieved.

The number of samples in the DOE was selected based on the cross-correlation between the variables. For each problem, the number of samples was increased until the cross-correlation between the variables became (\rho \le \mathrm{0.1}). In general, the number of samples is affected by several factors, including the number of design variables, the size of the FR, and the number of DOFs of the boundary variables. It appears that the latter factor has a greater impact on the resulting number of samples. For example, the first problem had four unknown parameters and two boundary variables with two DOFs (a translational spring at each end of the beam), and 500 samples were found to be adequate. The second example had six unknown parameters with four DOFs at the boundaries (a translational and a rotational spring at each end) that needed more samples; therefore, 1000 samples were used. The third problem had six unknown parameters but two unknown DOFs at the boundaries (two rotational springs); thus, 500 samples were deemed sufficient.

The order of the polynomials for the RSM depends on the magnitude of the coefficient of determination R², as specified in Eq. (2). For each problem, the order of the polynomial is increased until R² becomes≥0.99. It appears that the order of the polynomial is also sensitive to the number of DOFs at the boundaries. In the third example, only rotational DOFs were considered as variables; the second example was more sophisticated because both translational and rotational DOFs were considered in RVs. Therefore, the higher order of polynomials in the second example was required to reduce DOF interaction effects during RS construction and curve fitting.

The minimum value of FGP was selected by increasing the number of myriad sets (N_D). In this case, the feasible domain was discretized by a finer mesh and much smaller distances that provided a smaller FGP. Consequently, a smaller FGP returns system parameters with higher precision in their values. In addition, the FGP helps the experimenter become aware of the precision with which modal information (NFs and transmissibility) must be extracted in order to find the system’s unknown parameter with greater/lesser precision. In the first example, both (1%) and (5%) FGP were exploited to show how FGP affects the precision of the optimum points, as shown in Fig. 6. In the second and third examples, (0.1%) FGP was provided based on fine discretization of the FR by one million samples (myriad set) and the associated consecutive sample distances. As a rule of thumb, more samples produce higher resolution and result in higher precision of the optimum points.

The number of samples needed to obtain a myriad discrete set depends on the computer/LHS capacity and speed in generating samples as well as on the required accuracy of both system parameters and features. The more points used, the better, as they discretize the FR with a higher accuracy. Because of the limitations of our computing system, we were not able to go beyond one million points. Nevertheless, one million points produced (0.1%) FGP, which met the acceptable resolutions for the three examples. For instance, for the second NF of the second example and for the 5 Hz samples, [5.005 Hz] and [4.995 Hz] samples are available in the neighborhood. If analysts need higher resolutions of features, they can increase the number of myriad sets (N_D).

For the first example, nodes at \xi =\mathrm{0.2} and \xi =\mathrm{0.8} were chosen for the local information collection (MT) based on Eq. (6). These nodes attain relatively high STD while they are located away from the nodal points of the first three modes simultaneously, and they are also round numbers, which may make it easy for the analyst to use. Neighboring nodes at (\xi =\mathrm{0.18}, \xi =\mathrm{0.82}) and (\xi =\mathrm{0.22}, \xi =\mathrm{0.78}) were also selected, which represented deviations of up to 10% around \xi =\mathrm{0.2} and \xi =\mathrm{0.8}, but that did not affect the final values of the resulting system parameters. A similar approach was used to select the nodes for the second example at \xi =\mathrm{0.1} and at \xi =\mathrm{0.9}. In this case, neighboring nodes with (\xi =\mathrm{0.12}, \xi =\mathrm{0.88}) and (\xi =\mathrm{0.08}, \xi =\mathrm{0.92}) were also selected and did not affect the final values of the resulting system parameters.

There are still several situations in which the algorithm could struggle to find optimum points between the myriad set. First, either the PDE or FE model could be faulty and fail to recognize some system differential behaviors. Second, the EVP solver might return imprecise results, which would yield wrong RSs. Third, the RSs might not attain acceptable accuracy in estimating genuine structural responses (R^{\mathrm{2}}\le \mathrm{0.99}). Fourth, the myriad set might not be vast enough to contain genuine optimum points. Fifth, the FR might be not broad enough, and the system parameters might fall inside the pre-defined FR. Sixth, and most importantly, the experimentally extracted features may be fallacious, which can occur due to environmental setup and/or analyst error. Finally, selecting the locations of the critical node can be a challenging task, especially when a large number of mode shapes is used. In addition, these locations can be impractical for experimental implementation. Nevertheless, numerical experiments have shown that neighboring nodes are effective as long as they are far enough away from the nodal points of the Nm modes. To put it simply, if simultaneous high STDs of Nm local information can be achieved at neighboring nodes as well, those nodes can still be good candidates for data collection.

References

[1] Zong Z, Lin X, Niu J. Finite element model validation of bridge based on structural health monitoring—Part I: Response surface-based finite element model updating. Journal of Traffic and Transportation Engineering, 2015, 2(4): 258–278

[2] Castro-Triguero R, Murugan S, Gallego R, Friswell M I. Robustness of optimal sensor placement under parametric uncertainty. Mechanical Systems and Signal Processing, 2013, 41(1–2): 268–287

[3] Sehgal S, Kumar H. Structural dynamic model updating techniques: A state-of-the-art review. Archives of Computational Methods in Engineering, 2016, 23(3): 515–533

[4] Fang S E, Zhang Q H, Ren W X. An interval model updating strategy using interval response surface models. Mechanical Systems and Signal Processing, 2015, 60–61: 909–927

[5] Ren W X, Chen H B. Finite element model updating in structural dynamics by using the response surface method. Engineering Structures, 2010, 32(8): 2455–2465

[6] Zapico J L, González M P, Friswell M I, Taylor C A, Crewe A J. Finite element model updating of a small scale bridge. Journal of Sound and Vibration, 2003, 268(5): 993–1012

[7] Park Y S, Kim S, Kim N, Lee J J. Finite element model updating considering boundary conditions using neural networks. Engineering Structures, 2017, 150: 511–519

[8] Park Y S, Kim S, Kim N, Lee J J. Evaluation of bridge support condition using bridge responses. Structural Health Monitoring, 2019, 18(3): 767–777

[9] Cui Y, Lu W, Teng J. Updating of structural multi-scale monitoring model based on multi-objective optimisation. Advances in Structural Engineering, 2019, 22(5): 1073–1088

[10] Gordis J H. Artificial boundary conditions for model updating and damage detection. Mechanical Systems and Signal Processing, 1999, 13(3): 437–448

[11] Zhou L, Wang L, Chen L, Ou J. Structural finite element model updating by using response surfaces and radial basis functions. Advances in Structural Engineering, 2016, 19(9): 1446–1462

[12] Mottershead J E, Link M, Friswell M I. The sensitivity method in finite element model updating: A tutorial. Mechanical Systems and Signal Processing, 2011, 25(7): 2275–2296

[13] Jin S S, Cho S, Jung H J, Lee J J, Yun C B. A new multi-objective approach to finite element model updating. Journal of Sound and Vibration, 2014, 333(11): 2323–2338

[14] Jang J, Smyth A W. Model updating of a full-scale FE model with nonlinear constraint equations and sensitivity-based cluster analysis for updating parameters. Mechanical Systems and Signal Processing, 2017, 83: 337–355

[15] Nanthakumar S S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[16] Mao X, Dai H. A quadratic inverse eigenvalue problem in damped structural model updating. Applied Mathematical Modelling, 2016, 40(13–14): 6412–6423

[17] Tsai S H, Ouyang H, Chang J Y. Inverse structural modifications of a geared rotor-bearing system for frequency assignment using measured receptances. Mechanical Systems & Signal Processing, 2018, 110(Sep): 59–72

[18] Anitescu C, Atroshchenko E, Alajlan N, Rabczuk T. Artificial neural network methods for the solution of second order boundary value problems. Computers, Materials, & Continua, 2019, 59(1): 345–359

[19] Nehete D V, Modak S V, Gupta K. Experimental studies in finite element model updating of vibro-acoustic cavities using coupled modal data and FRFs. Applied Acoustics, 2019, 150: 113–123

[20] Catbas F N, Ciloglu S K, Hasancebi O, Grimmelsman K, Aktan A E. Limitations in structural identification of large constructed structures. Journal of Structural Engineering, 2007, 133(8): 1051–1066

[21] Zhang J, Maes K, De Roeck G, Reynders E, Papadimitriou C, Lombaert G. Optimal sensor placement for multi-setup modal analysis of structures. Journal of Sound and Vibration, 2017, 401: 214–232

[22] Mallardo V, Aliabadi M. Optimal sensor placement for structural, damage and impact identification: A review. Structural Durability and Health Monitoring, 2013, 9(4): 287–323

[23] Guo H Y, Zhang L, Zhang L L, Zhou J X. Optimal placement of sensors for structural health monitoring using improved genetic algorithms. Smart Materials and Structures, 2004, 13(3): 528–534

[24] Sun H, Büyüköztürk B. Optimal sensor placement in structural health monitoring using discrete optimization. Smart Materials and Structures, 2015, 24(12): 125034

[25] Jaishi B, Ren W X. Structural finite element model updating using ambient vibration test results. Journal of Structural Engineering, 2005, 131(4): 617–628

[26] Avitabile P. Modal Testing (A Practitioner’s Guide). Hoboken, NJ: John Wiley & Sons, 2018

[27] Marler R T, Arora J S. The weighted sum method for multi-objective optimization: New insights. Structural and Multidisciplinary Optimization, 2010, 41(6): 853–862

[28] Shahidi S G, Pakzad S N. Generalized response surface model updating using time domain data. Journal of Structural Engineering, 2014, 140(8): A4014001

[29] Rennen G, Husslage B, Van Dam E R, Den Hertog D. Nested maximin Latin hypercube designs. Structural and Multidisciplinary Optimization, 2010, 41(3): 371–395

[30] Brincker R, Ventura C. Introduction to Operational Modal Analysis. Hoboken, NJ: John Wiley & Sons, 2015

[31] Weijtjens W, Lataire J, Devriendt C, Guillaume P. Dealing with periodical loads and harmonics in operational modal analysis using time-varying transmissibility functions. Mechanical Systems and Signal Processing, 2014, 49(1–2): 154–164

[32] Devriendt C, Guillaume P. The use of transmissibility measurements in output-only modal analysis. Mechanical Systems and Signal Processing, 2007, 21(7): 2689–2696

[33] Devriendt C, Guillaume P. Identification of modal parameters from transmissibility measurements. Journal of Sound and Vibration, 2008, 314(1–2): 343–356

[34] Yan W J, Ren W X. An enhanced power spectral density transmissibility (EPSDT) approach for operational modal analysis: Theoretical and experimental. Engineering Structures, 2015, 102: 108–119

[35] Yan W J, Zhao M Y, Sun Q, Ren W X. Transmissibility-based system identification for structural health Monitoring: Fundamentals, approaches, and applications. Mechanical Systems and Signal Processing, 2019, 117: 453–482

[36] Brandt A. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Hoboken, NJ: John Wiley & Sons, 2011

[37] Marwala T. Finite element model updating using response surface method. In: Proceedings of the 45th Collection of Technical Papers-AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. California: Palm Springs, 2004, 5165–5173

[38] Box G E P, Draper N R. Empirical Model Building and Response Surface. Hoboken, NJ: John Wiley & Sons, 1987

[39] Myers R H, Montgomery D C, Anderson-Cook C M. Response Surface Methodology: Process and Product Optimization Using Designed Experiments. Hoboken, NJ: John Wiley & Sons, 2016

[40] Khuri A I, Mukhopadhyay S. Response surface methodology. Wiley Interdisciplinary Reviews: Computational Statistics, 2010, 2(2): 128–149

[41] Raseman W J, Jacobson J, Kasprzyk J R. Parasol: An open source, interactive parallel coordinates library for multi-objective decision making. Environmental Modelling & Software, 2019, 116: 153–163

[42] Huang T H, Huang M L, Jin J S. Parallel rough set: Dimensionality reduction and feature discovery of multi-dimensional data in visualization. In: Lu B L, Zhang L, Kwok J, eds. Neural Information Processing. ICONIP 2011. Lecture Notes in Computer Science, 7063. Berlin, Heidelberg: Springer, 2011.

[43] Inselberg A. The plane with parallel coordinates. Visual Computer, 1985, 1: 69–91

[44] Papazafeiropoulos G, Muñiz-Calvente M, Martínez-Pañeda E. Abaqus2Matlab: A suitable tool for finite element post-processing. Advances in Engineering Software, 2017, 105: 9–16

[45] Bishop R E D, Johnson D B. The Mechanics of Vibration. Cambridge: Cambridge University Press, 1960.

[46] Digilov R M, Abramovich H. Flexural vibration test of a beam elastically restrained at one end: A new approach for Young’s modulus determination. Advances in Materials Science and Engineering, 2013, 2013: 329530

[47] Leissa A W, Qatu M S. Vibration of Continuous Systems. New York: McGraw Hill Professional, 2011

[48] Fahy F, Walker J. Advanced Applications in Acoustics, Noise and Vibration. London: Taylor & Francis, 2005

[49] Gonçalves P J P, Brennan M J, Peplow A, Tang B. Calculation of the natural frequencies and mode shapes of a Euler-Bernoulli beam which has any combination of linear boundary conditions. Journal of Vibration and Control, 2019, 25(18): 2473–2479

[50] Li Z, Tang D, Li W. Analysis of vibration frequency characteristic for elastic support beam. Advanced Materials Research, 2013, 671–674: 1324–1328

[51] Karnovsky I A, Lebed O I. Free Vibrations of Beams and Frames: Eigenvalues and Eigenfunctions. New York: McGraw Hill Professional, 2004

[52] Dawson C B, Cha P D. A sensitivity-based approach to solving the inverse eigenvalue problem for linear structures carrying lumped attachments. International Journal for Numerical Methods in Engineering, 2019, 120(5): 537–566

[53] Devriendt C, De Sitter G, Vanlanduit S, Guillaume P. Operational modal analysis in the presence of harmonic excitations by the use of transmissibility measurements. Mechanical Systems and Signal Processing, 2009, 23(3): 621–635

[54] Araújo I G, Laier J E. Operational modal analysis using SVD of power spectral density transmissibility matrices. Mechanical Systems and Signal Processing, 2014, 46(1): 129–145

[55] Montgomery D C. Design and Analysis of Experiments. Hoboken, NJ: John Wiley & Sons, 2017

[56] Sansen J, Richer G, Jourde T, Lalanne F, Auber D, Bourqui R. Visual exploration of large multidimensional data using parallel coordinates on big data infrastructure. Informatics (MDPI), 2017, 4(3): 21

RIGHTS & PERMISSIONS

Higher Education Press