Robust improvement of the finite-element-model updating of historical constructions via a new combinative computational algorithm

Advances in Engineering Software 190 (2024) 103598 Available online 6 February 2024 0965-9978/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft Research paper Robust improvement of the finite-element-model updating of historical constructions via a new combinative computational algorithm Javier Naranjo-Pérez a, Rubén Rodríguez-Romero b, Pablo Pachón c,∗, Víctor Compán c, Andrés Sáez b, Aleksandar Pavic d, Javier Fernando Jiménez-Alonsob aDepartment of Continuum Mechanics and Structures, E.T.S. Ingenieros de Caminos, Canales y Puertos. Universidad Politécnica de Madrid, Madrid, Spain bDepartment of Continuum Mechanics and Structural Analysis, Universidad de Sevilla, Seville, Spain cDepartment of Building Structures and Ground Engineering, Universidad de Sevilla, Seville, Spain dVibration Engineering Section, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, UK ARTICLE INFO Keywords: Finite element model updating Maximum likelihood method Bi-objective optimization Decision-making problems Historical building stone ABSTRACT Finite-element-models are usually employed to simulate the behaviour of historical constructions. However, despite the high complexity of these numerical models, there are always discrepancies between the actual behaviour of the structure and the numerical predictions obtained. In order to improve their performance, an updating process can be implemented. According to this process, the value of the most relevant physical parameters of the model is adjusted to better mimic the actual behaviour of the structure. For this purpose, the actual structural behaviour is usually characterized via its experimental modal properties (natural frequencies and associated vibration modes). For practical engineering applications, the maximum likelihood method is normally considered to cope with this problem, due to its easy implementation together with an understandable interpretation of the updating results. However, the complexity of these numerical models makes unfeasible the practical implementation of the process due to the simulation time required for its computation. In order to shed some light to this problem, a new combinative computational algorithm is proposed herein. Additionally, the performance of the proposal has been assessed successfully via two applications: (i) a validation example, the model updating of a laboratory footbridge, in which the practical implementation of the algorithm has been described in detail; and (ii) a case-study, the model updating of a complex historical construction, in which the main advantage of the proposal has been highlighted, a clear reduction of the simulation time required to solve the updating problem without compromising the accuracy of the solution obtained. 1. Introduction Finite element (FE) method is usually employed to simulate numerically the structural behaviour of historical constructions [1]. However, despite the great complexity of these models [2], there are always differences between the predictions provided by these numerical tools and the actual behaviour of the structure [3]. These differences are normally originated by some of the following sources of errors and uncertainty [4]: (i) the approximate solution of differential equations that governs the dynamic behaviour of a structure (numerical uncertainty); (ii) the imprecise value of the model input parameters (parameter uncertainty); (iii) the incomplete definitions of underlying physics due to assumptions and idealizations (bias error); and (iv) the variability in measurements (experimental uncertainty). In order to reduce the effect of these sources of uncertainty, the FE model updating method is normally considered [5]. According to ∗Corresponding author. E-mail address: [email protected] (P. Pachón). this method, some physical parameters of the FE model are modified in order to better mimic the real behaviour of the structure. In the case of complex historical construction, the updated FE model is usually employed both to better assess the static or dynamic response of the structure and to establish a damage detection application based on a structural health monitoring strategy. Different criteria can be considered for the classification of the FE model updating methods. Regarding the time available for its computation, two types of methods can be considered: (i) off-line methods; and (ii) on-line methods. Thus, under the off-line approach, the FE model of the structure is directly updated via the modification of the value of its most relevant physical parameters. A high computational time is normally required to solve the updating problem according to his approach. This method is usually used to improve the accuracy of FE model when they are employed to assess the structural behaviour of https://doi.org/10.1016/j.advengsoft.2024.103598 Received 29 August 2023; Received in revised form 1 February 2024; Accepted 2 February 2024 Advances in Engineering Software 190 (2024) 103598 2 J. Naranjo-Pérez et al. civil engineering structures and historical constructions [2]. In addition, under the on-line approach, the FE model is approximated via a subrogate model and subsequently this model is calibrated according to some updating method. The use of this subrogate model allows reducing the simulation time required to solve the updating problem. This method is commonly used to damage detection under structural health monitoring strategies [6]. In this manner, the main difference between this two type of updating methods is related to the considered numerical model. Therefore, the same process is commonly employed to solve the updating problem for this two methods. This study is focused on the off-line updating approach, but its results can be easily extrapolated to the on-line approach. Additionally, FE model updating problems can be classified, regarding the parameter selection, in two main groups [7]: (i) direct methods (non-iterative) and (ii) indirect methods (iterative). On the one hand, according to the direct methods, the updating parameters are directly the components of the matrices that constitute the system equations considered to numerically mimic the behaviour of the structures. The main advantage of these methods is its direct and easy implementation but they present a clear drawback, the nonphysical meaning of the updated value of the parameters. On the other hand, according to the indirect methods, the updating parameters are some of the most relevant physical parameters of the structure (as for instance, Young’s modulus of the materials, stiffness of the support ...) that are iteratively modified in order to adjust the predictions of the numerical model to the actual response of the structure. In contrast to the previous methods, the adjusting process is more complicated, being necessary the use of complex computational tools, but they presents as advantage an easy physical interpretation of the updated value of the considered parameters. This advantage allows that the updating process can be considered a system identification method and it has originated that iterative methods have been widely used for the FE model updating of building and civil engineering structures. Among the indirect methods, a new classification can be established in terms of the quantification of the level of uncertainty. Two types of methods can be considered [7]: (i) deterministic and (ii) stochastics methods. Although, stochastics methods, as probabilistic [8] and fuzzy approaches, allow the determination of the level of uncertainty of the estimation of the updated values of the physical parameters [9], the high simulation time required to solve the updating problem, according to these approaches, have reduced, to date, their extensive practical implementation. For robust FE models of complex historical construction, deterministic methods are normally used. Among the deterministic methods, the maximum likelihood method has been widely employed due to its easy implementation. According to this method, the updating problem can be formulated as either a sensitivity-based weighted single-objective optimization problem or the combination of a bi-objective optimization sub-problem and a decisionmaking sub-problem. Some recent results, [10] about the performance of both methods when they are implemented for the FE model updating of civil engineering structures; have concluded that the second alternative allows obtaining a better adjustment since this second method explores in detail all the domain space of the possible value of the considered physical parameters. Therefore, the updating problem, according to this second method, is formulated considering a bi-objective function whose components are defined in term of the relative differences (residuals) between the experimental (actual) and numerical behaviour of the structure [11]. In order to characterize the actual behaviour of the structure, its experimental modal properties, identified considering either an experimental (EMA) [12] or operational modal analysis (OMA), are usually considered [13]. In this manner, these residuals are usually defined in terms of the relative differences between the experimental and numerical modal properties (natural frequencies and associated vibration modes) of the structure. Additionally, as the relation between the updating parameters and this bi-objective function is clearly nonlinear and, as a consequence, multiple local minimums can be reached; global optimization algorithms must be used to solve this problem [14]. As result of this optimization process, a set of possible solutions is obtained, the socalled Pareto front [12]. Finally, a subsequent decision-making problem must be solved, the selection of the best balanced solution (the ‘‘knee’’ point) among the different elements of the Pareto front. For this purpose, different methods have been considered [15] without existing to date a unified criterion to cope with this problem. In spite of all the advantages of the abovementioned method, when it is implemented to solve the FE model updating of complex historical construction structures, it presents a clear drawback, the high simulation time required to solve the problem. In order to improve the performance of this method, several numerical techniques can be adopted: (i) to increase the search speed of the optimization algorithm combining local and global computational optimization algorithms [16]; (ii) to reduce the complexity of the FE model approximating its behaviour via surrogate models [17]; (iii) to find the optimum selection of the physical parameters of the FE model; (iv) to use alternative mathematical tools to solve the updating problem more efficiently [7]; and (v) to use some clustering technique for the parallel computation of the updating problem [18]. This research study focuses in the first alternative, the proposal of a new combinative computational algorithm to reduce the simulation time required to perform the FE model updating of historical construction without compromising the accuracy of the solution obtained. For this purpose, this proposal takes advantage of the virtues of two previous algorithms [19,20], proposed by some of the authors of the manuscript, and it combines them to improve the efficiency of the updating process. For this purpose, the proposed algorithm makes use of a local– global optimization algorithm, the UKF-MHS algorithm [19], both to increase the convergence speed of the searching algorithm and to reduce the uncertainty associated with the variability of the experimental measurements, together with a collaborative algorithm [20], which allows a robust estimation of the ‘‘knee’’ point via the linking of several statistical learning techniques. The combination between these two algorithms improves clearly the performance of the FE model updating of complex historical constructions. Thus, it is necessary to highlight that the main contribution of this paper is not only the combination of two previously reported algorithms to obtain an efficient computational tool that assists architects and structural engineers in the FE model updating of historical constructions but also its practical implementation in two examples. These application examples try to highlight the virtues of the proposal. First, this new algorithm has been implemented for the off-line FE model updating of a laboratory footbridge. In this manner, this first application example has been used as tutorial to explain in detail the different steps of this computational algorithm. Subsequently, the algorithm has been implemented for the off-line FE model updating of a complex historical constructive. In this second application example, the higher performance of the proposal in comparison with the abovementioned conventional bi-objective approach has been remarked. In addition to this introductory section, the remaining paper is organized as follows. Section 2describes both the basics of the FE model updating, according to the maximum likelihood method, and the formulation of the proposed new combinative algorithm. Section 3 focuses on the detailed description of the proposal via a validation example, the practical implementation of the algorithm for the FE model updating of a laboratory footbridge. Section 4presents the performance assessment of this proposal when it is implemented for the FE model updating of a historical construction. For this purpose, the model updating of the church of the Royal Monastery of San Jerónimo (Granada, Spain) is developed in detail. Finally, Section 5summarizes the main conclusions that can be drawn from the work. Advances in Engineering Software 190 (2024) 103598 3 J. Naranjo-Pérez et al. 2. Solving the maximum likelihood FE model updating problem via a new combinative computational algorithm 2.1. Basics of FE model updating according to the maximum likelihood method The goal of FE model updating is to acquire a precise numerical FE model that corresponds to the experimental dynamic properties of the analysed civil or building engineering structure. This objective can be accomplished by approximating the most significant physical parameters of the model that reduce the differences between the numerical and experimental data (residuals). Consequently, the FE model updating issue can be interpreted as a parameter identification problem. In this regard, the maximum likelihood method is the most commonly utilized estimator for civil or building engineering applications. Assuming that the residuals follow a normal distribution, the maximum likelihood method is analogous to the ordinary least squares method. Hence, the model updating can be converted into an optimization problem, whose primary goal is to reduce the total of the relative deviations between the experimental and numerical modal properties. Thus, the updating problem can be formulated alternatively as two different optimization problems [7]: (i) a sensitivity-based weighted single-objective optimization problem; or (ii) a combination of a bi-objective optimization sub-problem together with a decision-making sub-problem. Due to the higher efficiency of the second formulation [10] to cope with this problem, the multi-objective approach has been considered herein. The relevant physical parameters of the model are taken as the design variables for this purpose. Therefore, the formulation of the FE model updating issue considering the bi-objective approach can be represented as: min [𝑓1(𝜽)𝑓2(𝜽)] = min [ 1 2 𝑚𝑓 ∑ 𝑗=1 𝑟𝑓 𝑗(𝜽)21 2 𝑚𝑚 ∑ 𝑗=1 𝑟𝑚 𝑗(𝜽)2] 𝜽𝑙≤𝜽≤𝜽𝑢 (1) where 𝑟𝑓 𝑗(𝜽) and rmj 𝑟𝑚 𝑗(𝜽) refer to the discrepancies between the 𝑗th natural frequency and mode shapes of the structure and their theoretical counterparts; 𝑚𝑓represents the quantity of considered natural frequencies; 𝑚𝑚represents the quantity of considered mode shape; 𝜽is a vector comprising the updating parameters of the FE model; and 𝜽𝑙 and 𝜽𝑢denote the lower and upper limits of the exploration range for these physical parameters, respectively. The residuals may be defined as: 𝑟𝑓 𝑗(𝜽) = 𝑓𝑛𝑢𝑚,𝑗 (𝜽) − 𝑓𝑒𝑥𝑝,𝑗 𝑓𝑒𝑥𝑝,𝑗 𝑗= 1,2,…, 𝑚𝑓(2) 𝑟𝑚 𝑗(𝜽) = √ √ √ √ √⎛ ⎜ ⎜ ⎝ (1 − √𝑀𝐴𝐶,𝑗 (𝜽))2 𝑀𝐴𝐶,𝑗 (𝜽)⎞ ⎟ ⎟ ⎠ 𝑗= 1,2,…, 𝑚𝑚(3) considering the numerical natural frequency 𝑗as 𝑓𝑛𝑢𝑚,𝑗 (𝜽), the experimental natural frequency 𝑗as 𝑓𝑒𝑥𝑝,𝑗 (𝜽), and the Modal Assurance Criterion value 𝑀𝐴𝐶,𝑗 (𝜽)[21] as a measure used to assess the similarity between the mode shapes of the numerical vibration mode 𝑗,𝜑𝑛𝑢𝑚,𝑗 (𝜽) and the experimental vibration mode 𝑗,𝜑𝑒𝑥𝑝,𝑗 (𝜽). As the accuracy of the 𝑀𝐴𝐶,𝑗 (𝜽)ratio depends on the geometrical definition of the mode shapes, a sufficient refined gridline is needed to characterize the value of this magnitude [10]. For this purpose, a sensitivity analysis can be performed. In this study, the variation of the 𝑀𝐴𝐶,𝑗 (𝜽)ratio in term of the number of coordinates that defined the gridline (of the experimental identification test) is analysed in detail. A balanced value between the complexity of the gridline and the accuracy of the 𝑀𝐴𝐶,𝑗 (𝜽)ratio is taken into account. The physical parameters of the updating process must be selected carefully since they condition the accuracy and reliability of the solution. Based on a FE model, a different updating problem can be formulated depending on the relation between the residuals and the considered updating parameters. For the parameter selection purpose, different criteria can be considered [7]. Among the different criteria, a sensitivity analysis is commonly employed. According to this study, the variation of the natural frequencies and mode shapes in term of the considered physical parameters is computed. As this variation is proportional to the modal strain energy this magnitude is considered to perform the sensitivity analysis. The physical parameters, which reflect a higher variation, are selected as updating parameters [10]. The resolution of this bi-objective optimization problem generates a set of parameter vectors, each representing a potential solution. These parameters can be visually represented by forming a curve known as the Pareto front. Each point along the Pareto front represents an optimal updated model where improving one objective would require sacrificing another. Consequently, selecting the best solution from the Pareto front requires a decision-making process. Several methods have been proposed in literature to tackle this challenge by balancing the different terms of the Pareto front. Generally, these methods identify the best solution on the Pareto front as the point where a slight improvement in one objective would significantly deteriorate the other objective. This particular point on the Pareto front is often referred to as the ‘‘knee’’ point, and each method suggests a slightly different criterion for its determination. Following the above approach and under the maximum likelihood method, the FE model updating is usually addressed by implementing computational intelligence algorithms [22]. For this purpose, a FE analysis package is commonly linked to an optimization package. As illustration of this connectivity, Fig. 1 shows the flowchart of the updating process giving a special emphasis to the linking between the two abovementioned packages together with the flow of data. During the updating process, a pairing problem can be solved. The numerical and experimental mode shapes must be matched. For this purpose, the 𝑀𝐴𝐶𝑗(𝜽)ratio is commonly considered. For this purpose, the 𝑀𝐴𝐶𝑗(𝜽)ratio among all the numerical and experimental mode shape is computed and organized in a matrix. The rows of this matrix are associated with the numerical mode shapes and the columns with the experimental mode shapes. The component of the matrix with a higher value allows determining the linking among numerical and experimental vibration modes. Although computational algorithms have been used extensively to cope with this updating problem, they present two clear limitations. First, they elapse a high simulation time when the complexity of the FE model increases. As it is the case of robust FE models. Second, they are not able to deal with uncertainties associated with the experimental modal properties of the system. To tackle these two issues, a hybrid local–global algorithm combined with a predictive model is used in this paper. In particular, the hybrid Unscented Kalman Filter (UKF)-Harmony Search (HS) algorithm is used to obtain the Pareto front of the bi-optimization problem [19, 23]. This algorithm allows considering the uncertainties of both the measurements and the estimation process and takes advantage of the acceleration scheme provided by the UKF and the global character of the HS. Once the Pareto front is obtained, a decision-making problem must be solved to select the best solution. To cope with this, the Pareto front is processed through a PC analysis and an artificial neural network (ANN) is designed to map the relationship between the physical parameters of the model and the bi-objective function. Finally, the best solution is determined from a local minimization problem. In this sense, the ANN defines a continuous predictive model that allows reducing the number of iterations, thus, the simulation time of the optimization algorithm [20]. In this section, the multi-objective HS, the UKF algorithm and the post-process of the Pareto front are detailed. Advances in Engineering Software 190 (2024) 103598 4 J. Naranjo-Pérez et al. Fig. 1. Flowchart representing the linking between the FE analysis package and the mathematical software for the FE model updating implementation. 2.2. MHS algorithm The global metaheuristic HS algorithm [24] is inspired in the mental process carried out by musicians when they improvise searching for harmony according to aesthetic requirements. The goal of this algorithm is to find the global minimum of an objective function via the iterative modification of a set of physical parameters of the model. When compared to Genetic Algorithm, the complexity of the mathematical operations implemented in the HS algorithm to monitor the evolution of the population are designed to improve its performance. The HS algorithm has been implemented successfully for several practical engineering applications [25]. This algorithm has shown a great effectiveness when compared with other conventional metaheuristic algorithms to solve the FE model updating problem of complex civil engineering structures [10]. The HS algorithm consists of the following steps [10]. First, the harmony matrix, H, that stores the initial candidate solutions is created. After that, for each candidate solution the objective function is assessed. The new harmonies are generated using three mechanisms (memory consideration, pitch adjustment and randomization) and the objective function is evaluated again. Finally, the harmony matrix, H, is updated with the best harmonies and the steps are repeated until a convergence criterion is reached. The MHS is an extension of the above described algorithm which allows minimizing multi-objective functions. When a new harmony is generated, each element of a new candidate vector can be defined in terms of either a previous value stored in the harmony matrix, H, or adopting a new random value. This fact is controlled by the harmony memory consideration rate, HMCR, Thus, for each parameter,a random number, between 0 and 1, is generated and if this number is lower than or equal to HMCR, the parameter is chosen randomly from the matrix H. Otherwise, a random value among the possible values of the search domain is assigned to the parameter. Additionally, when some elements adopts the value of a previous one, it can be mutated according to the pitch adjustment rate, PAR. This parameter works in the same way as the previous one. In case the parameter must be adjusted, this is based on an additional parameter, the so-called bandwidth, b𝑤, which is added or subtracted to mutate the considered candidate vector. The classification of the non-dominated solutions is performed using both the non-dominated sorting method and the crowding distance [26]. Finally, in order to restore the initial size of the harmony matrix, H, the worst solutions, according to a crowding distance criterion, are removed. When the convergence criterion is met, the so-called Pareto front is obtained in terms of the non-dominated solutions. 2.3. Unscented Kalman filter The Unscented Kalman Filter (UKF) is framed within the family of sigma-points Kalman filters [27]. This local algorithm is a derivative free estimator, i.e., no Jacobians or Hessians need to be calculated, broadly employed for nonlinear dynamic systems. The formulation for a parameter estimation problem is represented as [28]: 𝜽𝑘=𝜽𝑘−1 +𝐰𝑘−1 (4) 𝐳𝑘=ℎ(𝜽𝑘) + 𝐯𝑘(5) being 𝜽the parameter vector; ℎthe nonlinear modelling function, 𝐳the outputs of the dynamic system; 𝐰the statistical noise of the identification process; and 𝐯the statistical noise of observation process. Both type of noise are assumed to be uncorrelated and white Gaussian noise with zero-mean and covariance matrices 𝑸and 𝑹respectively. The matrix 𝑹may be computed by means of two terms [29]: (i) the measurement noise; and (ii) the modelling noise. If the same FE model is considered for each step, the first component of the matrix 𝑹can be neglected [16]. This algorithm performs the nonlinear estimation through the definition of 2𝑛+ 1 (being 𝑛the number of parameters) deterministic sampling points (sigma points). Their propagation through the nonlinear function, ℎ, leads to the true posterior mean and covariance of the estimated parameters up to the second order of the Taylor series expansion of the nonlinear function (third order of the Taylor series expansion for a Gaussian inputs). Hence, it is an extension of the unscented transformation [30]. The main computational effort of this algorithm is the computation of the new sigma points which is based on the square-root decomposition of the posterior covariance matrix, 𝐏. Thus, this matrix must be positive semidefinite at each step. To employ the Cholesky factorization (𝐀=√𝐏=𝑐ℎ𝑜𝑙(𝐏)being 𝐏=𝐀𝐀𝑇) may provide an efficient decomposition but the matrix 𝐏needs to be updated at each iteration and numerical errors can give a nonpositive semidefinite matrix. To overcome this issue, the square-root UKF [31] propagates directly the matrix 𝐴avoiding the factorization and ensuring that the matrix is positive semi-definite. These sigma points and their associated weights are calculated as: (𝝌𝑘−1)0= 𝜽𝑘−1|𝑘−1 (6) (𝝌𝑘−1)𝑖= 𝜽𝑘−1|𝑘−1 +√(𝑛+𝜆)(𝐀𝜃 𝑘−1|𝑘−1)𝑖𝑖= 1,2,…, 𝑛 (7) Advances in Engineering Software 190 (2024) 103598 5 J. Naranjo-Pérez et al. (𝝌𝑘−1)𝑖+𝑛= 𝜽𝑘−1|𝑘−1 −√(𝑛+𝜆)(𝐀𝜃 𝑘−1|𝑘−1)𝑖𝑖= 1,2,…, 𝑛 (8) 𝑊0=𝜆 𝑛+𝜆(9) 𝑊𝑖=𝑊𝑖+𝑛=1 2(𝑛𝑑+𝜆)𝑖= 1,2,…, 𝑛 (10) where  𝜃𝑘−1|𝑘−1 are the posterior parameters estimated at the previous step and 𝜆is a scaling parameter. As an algorithm belonging to the Kalman Filter family, the squareroot UKF has also two steps: prediction and correction. The former considers the prior model to assess the sigma points, predict the estimation error covariance, 𝐀𝜃, the model outputs error covariance, 𝐒𝑧, and the cross covariance between the estimation error and the model outputs error covariances, 𝐏𝜃𝑧. The estimation error covariance is given by 𝐀𝜃 𝑘=𝛾−1∕2𝐀𝜃 𝑘−1, where 𝛾is a scalar factor slightly less than the unit [31]. Since 𝑊𝑖>0for all 𝑖≥1, the model output error covariance, 𝐒z, can be expressed as [32]: 𝐒z 𝑘= 2𝑛𝑑 ∑ 0 𝑊𝑖[[(𝑧𝑖 𝑘|𝑘−1 − z𝑘|𝑘−1)] ⋅(𝑧𝑖 𝑘|𝑘−1 − z𝑘|𝑘−1)𝑇] + 𝐑 = [√𝑊𝑖(𝑧𝑖 𝑘|𝑘−1 − z𝑘|𝑘−1),√𝐑] ⋅[√𝑊𝑖(𝑧𝑖 𝑘|𝑘−1 − z𝑘|𝑘−1)𝑇,√𝐑 𝑇 ]𝑇 +𝑊0[(𝑧0 𝑘|𝑘−1 − z𝑘|𝑘−1)⋅(𝑧0 𝑘|𝑘−1 − z𝑘|𝑘−1)𝑇]𝑓𝑜𝑟 𝑖 =1∶𝑛𝑑 (11) The first term can be expressed by means of a QR factorization. The last term, can be taken into account performing a rank 1 update to Cholesky factorization. The cross covariance is then calculated as: 𝐏𝜃z 𝑘|𝑘−1 = 2𝑛𝑑 ∑ 0(𝑊𝑖[(𝝌𝑘−1)𝑖− 𝜽𝑘|𝑘−1]⋅[(𝑧𝑘|𝑘−1)𝑖− z𝑘|𝑘−1]𝑇)(12) Finally, the correction step uses both the measurements, z𝑜𝑏𝑠, and the model outputs to determine the posterior mean and covariance of the parameter estimation by considering the Kalman’’s gain matrix, 𝐊. 𝐊𝑘= (𝐏𝜃z 𝑘|𝑘−1∕𝐒z 𝑘|𝑘−1 𝑇)∕𝐒z 𝑘|𝑘−1 (13)  𝜽𝑘|𝑘= 𝜽𝑘|𝑘−1 +𝐊𝑘(z𝑜𝑏𝑠 − z𝑘|𝑘−1)(14) 2.4. Post-process of the pareto front The Pareto front is projected from the functional space to the principal component (PC) space by performing a principal component analysis (PCA) of the constituent points. PCA is a statistical tool to extract patterns from correlated variables following the next steps [33]: first, the data are normalized and their covariance matrix is computed; subsequently, the eigenvalues and eigenvectors of this matrix are calculated, and finally, the PC scores are determined. It is worth pointing out that the eigenvector associated to the highest eigenvalue is the first PC and so on and that the scores represent the projection of the data in the PC space. The resulting Pareto front exhibits a clear convex character which eases the determination of the best solution among all the possible ones as this is addressed solving a local minimization problem. The approximation of the previous Pareto front to a continuous predictive model is achieved employing an artificial neural network (ANN). This is a supervised machine learning technique that maps nonlinear and coupled relationships between a set of inputs and outputs. The multilayer perceptron (MLP) is the most widely used network topology and consists of an input layer, an output layer and one or several hidden layers, which are interconnected in a feed-forward manner: one neuron sends information to the subsequent layer but does not receive from them. The number of hidden layer can be assumed to be equal to one as it has been demonstrated that only a single hidden layer may approximate any continuous function [34,35]. A trial and error procedure or empirical relationships may me employed to determine the number of neurons. The application of the ANN to the Pareto front allows mapping the relationship between the physical parameters of the model and the multi-objective function playing the role of a convex function. The relationship needs to be trained and tested on the basis of a backpropagation training algorithm that adjusts the weights among the neurons’ connection by minimizing the error between the estimated output and the desired output propagated in a backward direction [36]. Taking advantage of the convex character of the continuous function, the decision-making problem to obtain the best solution is transformed into a local minimization problem due to its convergence speed and good accuracy. The solution obtained is the knee-point of the original Pareto front. As local optimization algorithm, the Active-Set (A-S) is hereby adopted [37]. Under this approach, the decision-making problem may be formulated as: minimize𝑔(𝜃) = 𝑦(𝜃)(15) Subjectedto ∶ ⎧ ⎪ ⎨ ⎪ ⎩ 𝜃𝑙≤𝜃≤𝜃𝑢 𝑥(𝜃)≥𝑥𝑚𝑖𝑛 𝑦(𝜃)≥𝑦𝑚𝑖𝑛 (16) being 𝑔(𝜃)the objective function, 𝜃𝑙and 𝜃𝑢are respectively the minimum and maximum values of the physical parameters which form the Pareto front, (𝑥(𝜃), 𝑦(𝜃)) are the ANN outputs in the PC space, and 𝑥𝑚𝑖𝑛 and 𝑦𝑚𝑖𝑛 are the minimum values of the ANN outputs in the PC space. 2.5. Proposed algorithm The implemented algorithm combines the main virtues of the two described algorithms (MHS and UKF) and the techniques employed to post-process the Pareto front. First, the uncertainties are considered and, second, the computational cost is reduced due two reasons. The first one is the acceleration scheme of the hybrid algorithm itself. The second one lies in the ANN which allows reducing the iterations and the population to obtain a non-crowded Pareto front as this is then simulated by the ANN. Finally, the PCA improves the efficiency of the decision-making problem since it provides a convex representation of the Pareto front which is further employed by the A-S algorithm to obtain the best solution in a straightforward manner (Fig. 2). 3. Validation example: FE model updating of a laboratory footbridge A real laboratory footbridge is considered as benchmark structure to validate the proposed hybrid-collaborative algorithm. Section 3.1 describes the geometry and main constitutive elements of the structure. Section 3.2 introduces the initial FE model of the footbridge. Then, Section 3.3 presents the dynamic characterization of this structure via a forced vibration test (FVT) together with an EMA algorithm. Subsequently, Section 3.4 details the parameter selection of for the updating process. Finally, Section 3.5 describes in detail the updating process under the maximum likelihood method. 3.1. Geometrical configuration This benchmark structure is a reconfigurable steel footbridge from the Vibration Engineering Section of the University of Exeter (Fig. 3). The footbridge is a frame structure with a length of 15 m. The structure is comprised of two lateral steel beams separated transversally 2.5 m. Rectangular plates of 200 ×12 mm spaced 1.25 m brace these two beams. Steel bolts connect the deck, made of composite SPS panels [38], with both the longitudinal beams and transversal plates. Four columns directly pinned to the ground are used as support of the four ends of the steel beams. Further information about this footbridge can be found in Ref. [39]. Advances in Engineering Software 190 (2024) 103598 6 J. Naranjo-Pérez et al. Fig. 2. Flowchart of the proposed combinative computational algorithm. Fig. 3. Description of the laboratory footbridge and layout of the FVT. 3.2. Initial FE model An initial FE model of the footbridge was built using the Ansys Software [40]. To model the structure, three different types of elements were employed: (i) shell elements (SHELL181) for the lateral beams; the transversal beams and the SPS panel; (ii) 3D beam elements (BEAM188) for the bolts of the connection between the SPS panel and the steel structure; and (iii) equivalent spring element (COMBIN14) with longitudinal and transversal stiffness to model each support whose vertical displacement was constrained. The value of the stiffness at each direction was estimated from a simplified FE model of just the column, resulting in an equivalent stiffness of 5.5⋅107N∕m and 1.0⋅107N∕m for the longitudinal and transversal direction, respectively. The mesh consisted of 31 903 elements. The adopted mechanical properties for Advances in Engineering Software 190 (2024) 103598 7 J. Naranjo-Pérez et al. Fig. 4. Sensitivity study to select the most relevant physical parameters of the footbridge model. Table 1 Experimental natural frequency, 𝑓𝑒𝑥𝑝 , initial numerical natural frequency 𝑓𝑖𝑛𝑖, relative difference, ▵𝑓, and the MAC ratio. Vibration mode 𝑓exp [Hz] 𝑓ini [Hz] ▵𝑓[%] MAC [-] 1 3.810 3.638 −4.51 0.998 2 5.144 5.433 5.62 0.994 3 8.485 9.106 7.32 0.988 4 12.366 11.310 −8.54 0.877 5 18.605 17.364 −6.67 0.985 6 20.459 19.519 −4.59 0.993 7 22.980 20.725 −9.81 0.634 the different materials are the followings: (i) for the steel, a density of𝜌𝑠= 7850 kg∕m3, the Young’s modulus equals 𝐸𝑠= 210 ⋅109Pa and a Poisson’s ratio of 𝑣𝑠= 0.3; and (ii) for the polyurethane the density was assumed to be 𝜌𝑝= 1100 kg∕m3, the Young’s modulus 𝐸𝑝= 7.5⋅108 Pa and the Poisson’s ratio 𝑣𝑝= 0.3. The numerical modal analysis of this model gives as result the seven natural frequencies shown in Table 1 and the mode shapes depicted in Fig. 5. 3.3. FVT and EMA The experimental modal identification of the footbridge was carried out through the experimental modal analysis of the accelerations recorded in a forced vibration test. To do this, a set of proof-mass actuators and accelerometers were employed (Fig. 3). Random signals were used as input to simultaneously drive the actuators within a Multiple Input-Multiple Output (MIMO) configuration [41]. From the applied forces and accelerations response, the Frequency Response Function was calculated considering an overlap of 50%. After this, a complex mode indicator function was defined to identify probable mode location in the FRFs fitted curves. To finish the process, a global polynomial curve fitting method extracts the first seven experimental natural frequencies given in Table 1 and their mode shapes (shown in Fig. 5) from the abovementioned probable mode locations. The reader can refer to Ref. [39] for further explanation regarding the forced vibration test and the experimental modal analysis. 3.4. Parameters selection and search domain of the FE model updating process Based on the differences between the numerical and experimental modal properties, the FE model of the footbridge is updated applying the proposed algorithm. Firstly, the physical parameters of the model that are being updated must be identified. The ratio between the modal strain energy associated with each parameter and the overall modal strain energy is considered as criterion to reveal the influence of the parameter on the variation of the natural frequencies [42]. Thus, the parameters with greater ratios are selected. An initial set of 15 physical parameters were taken into account in this analysis and, as result, ten were selected as design parameters (designated as in Fig. 4): Young’s modulus of the steel of the longitudinal beams at 6 different sections (𝜃𝑖𝑛,1 − 𝜃𝑖𝑛,6), Young’s modulus of the polyurethane (𝜃𝑖𝑛,8), Young’s modulus of the steel of the bolts (𝜃𝑖𝑛,9), equivalent longitudinal stiffness of the supports (𝜃𝑖𝑛,11) and the equivalent transversal stiffness of the supports (𝜃𝑖𝑛,12). The selection of the Young’s modulus of the different material as physical parameter has not got as objective the identification of this structural property via the resolution of the inverse problem. However, this magnitude is used as estimator of the stiffness of the structure. In this manner, the variation of this quantity with respect to its reference value reflects the reduction of some uncertainties (geometrical tolerances, constitutive laws ...) associated with the modelling of the structure. In addition, a search domain has been established for each design parameter in order to constrain the optimization problem and to guarantee an adequate physical meaning of the updated value. The considered research domain for each considered updating parameter with its lower and upper bound is illustrated by Table 2. Once both the design parameters and their search range have been set, the FE model updating process is performed. 3.5. FE model updating process The previous described bi-objective algorithm is applied for the model updating of the footbridge. The followings values for the parameters of the MHS algorithm were adopted [22,43]: a HMCR ratio of 0.9, a PAR ratio of 0.7 and a bandwidth, 𝑏𝑤, equals the 1% of the search domain of each parameter. Regarding the parameters of the UKF algorithm, the sensitivity analysis carried out by the authors in Ref. [19] was used, thus, the values of the parameters were: number of iterations of the UKF algorithm, NUKF=3, initial estimation error covariance, 𝑃𝜃 0=𝑑𝑖𝑎𝑔(((𝜃𝑢−𝜃𝑙)∕2000)2)and measurement noise covariance matrix, 𝑅𝑖𝑖 = 0.001. The number of iterations of the MHS algorithm was set to 10, the number of initial harmonies was 20 and the number of new harmonies, generated at each iteration, equals 5. As result, the non-crowded Pareto front containing the possible solutions was obtained (see Fig. 6a). The computational cost of this step of the Advances in Engineering Software 190 (2024) 103598 8 J. Naranjo-Pérez et al. Fig. 5. First seven mode shapes obtained from: (a) the experimental identification test; and (b) the initial FE model. proposed algorithm is the highest, as it took 4678 s to obtain the non-crowded Pareto front. Next, the PCA analysis of this Pareto front is conducted according to the steps of the proposed algorithm. The singular value decomposition method was employed to perform the decomposition of the covariance matrix. Fig. 6b shows the projection of the Pareto front in the PC space. As mentioned in Section 2, this analysis presents two advantages. Firstly, the accuracy of the ANN used in the next step to approximate the Pareto front is improved and, secondly, the selection of the best solution is more robust due to the convex nature of the Pareto front in the PC space. The time to perform this analysis is negligible when compared to the bi-objective optimization as it takes less than one second. The fourth step of the algorithm involved the design of the ANN to approximate the Pareto front. The objective is to simulate the relationship between the ten physical parameters and the two residuals of the bi-objective function in a continuous manner, yielding a continuous Pareto front which leads to a simpler decision-making problem. In this study the MLP (Multi-layer perceptron) with one hidden layer has been used as ANN topology. The rule of Kermanshahi [44] allows calculating the number of neurons of the hidden layer as neurons = (𝑚+𝑛)∕2 + 𝛿, being 𝑚and 𝑛the number of neurons of the input and output layers, respectively, and 𝛿a normalization factor which can be equal to 1 or 2. Thus, for this case 𝑚= 10,𝑛= 2 and 𝛿= 1, giving as result 7 neurons for the hidden layer. The ANN was trained using the Levenberg–Marquardt backpropagation algorithm [36], defining the error function in terms Advances in Engineering Software 190 (2024) 103598 9 J. Naranjo-Pérez et al. Fig. 6. Combinative computational algorithm: (a) non-crowded Pareto front, (b) PCA analysis, (c) solution using the ANN and (d) corresponding solution in the original Pareto front. Table 2 Updated value of the physical parameters of model, 𝜃𝑖𝑛, after the updating process. Parameter Description Initial value Range of variation Updated value Lower Upper 𝜃𝑖𝑛,1𝐸𝑠long. beam section 1 [GPa] 210 190 240 230.08 𝜃𝑖𝑛,2𝐸𝑠long. beam section 2 [GPa] 210 190 240 215.94 𝜃𝑖𝑛,3𝐸𝑠long. beam section 3 [GPa] 210 190 240 202.43 𝜃𝑖𝑛,4𝐸𝑠long. beam section 4 [GPa] 210 190 240 215.92 𝜃𝑖𝑛,5𝐸𝑠long. beam section 5 [GPa] 210 190 240 193.52 𝜃𝑖𝑛,6𝐸𝑠long. beam section 6 [GPa] 210 190 240 214.58 𝜃𝑖𝑛,8𝐸𝑝polyurethane [MPa] 1000 500 1500 758.92 𝜃𝑖𝑛,9𝐸𝑠steel bolts [GPa] 1000 210 2100 790.95 𝜃𝑖𝑛,11 KEquivalent long. stiffness [N/m2]6.0⋅1071.4⋅1071.1⋅1087.06 ⋅107 𝜃𝑖𝑛,12 KEquivalent transv. stiffness [N/m2]9.0⋅1064.8⋅1063.8⋅1077.79 ⋅106 of the mean square error and considering 70% of the elements of the Pareto front for sampling. The rest of the elements were used to validate and test the ANN. Fig. 6b illustrates the accuracy of the ANN used to continuously approximate the Pareto front in the PC space. Finally, the decision-making problem is solved to obtain the best solution (knee point). For this purpose, the A-S algorithm has been used to conduct the minimization problem. Once this point is obtained, it is projected back into the original Pareto front (Fig. 6c-d). The simulation time to design the ANN and calculate the knee point was about 5 s. The physical parameters associated to this optimum point (updated values) are collected in Table 2 and the corresponding modal properties of the updated FE model are shown in Table 3. It can be observed that after the updating process, all the MAC ratios are above 0.89 and the relative differences have been reduced. In addition, the computational cost to perform this process has been also less than traditional bi-objective optimization as it is not needed Table 3 Experimental natural frequency, 𝑓𝑒𝑥𝑝 , updated numerical natural frequency 𝑓𝑢𝑝𝑑 , relative difference, ▵𝑓, and the MAC ratio. Vibration mode 𝑓exp [Hz] 𝑓upd [Hz] 𝛥𝑓 [%] MAC [-] 1 3.810 3.844 0.892 0.999 2 5.144 5.458 6.104 0.994 3 8.485 8.388 1.143 0.988 4 12.366 11.858 4.107 0.890 5 18.605 18.148 2.456 0.986 6 20.459 20.185 1.339 0.993 7 22.980 21.607 5.975 0.963 a populated Pareto front, and, therefore, the number of iterations and population of the algorithm can be set to low values. Advances in Engineering Software 190 (2024) 103598 16 J. Naranjo-Pérez et al. [24] Geem ZW, Kim JH, Loganathan G. A new heuristic optimization algorithm: Harmony search. Simulation 2001;76(2):60–8. http://dx.doi.org/10.1177/ 003754970107600201. [25] Yang X-S, Koziel S, editors. Computational optimization and applications in engineering and industry. Springer Berlin Heidelberg; 2011, http://dx.doi.org/ 10.1007/978-3-642-20986-4. [26] Deb K, Agrawal S, Pratap A, Meyarivan T. A fast and elitist multiobjective genetic algorithm: Nsga-ii. IEEE Trans Evol Comput 2002;6:182–97. [27] Julier SJ, Uhlmann JK. New extension of the Kalman filter to nonlinear systems. In: Defense, security, and sensing. 1997. [28] Wan E, Van Der Merwe R. The unscented Kalman filter for nonlinear estimation. In: Proceedings of the IEEE 2000 adaptive systems for signal processing, communications, and control symposium(Cat. no. 00EX373). 2000, p. 153–8. http://dx.doi.org/10.1109/ASSPCC.2000.882463. [29] Tarantola A. Inverse problem theory and methods for model parameter estimation. Society for Industrial and Applied Mathematics; 2005, http://dx.doi.org/ 10.1137/1.9780898717921. [30] Julier S, Uhlmann J. Corrections to unscented filtering and nonlinear estimation. Proc IEEE 2004;92(12):1958. http://dx.doi.org/10.1109/JPROC.2004.837637. [31] Van der Merwe R, Wan E. The square-root unscented Kalman filter for state and parameter-estimation. In: 2001 IEEE international conference on acoustics, speech, and signal processing. proceedings (Cat. no. 01CH37221), vol. 6. 2001, p. 3461–4. http://dx.doi.org/10.1109/ICASSP.2001.940586. [32] Terejanu G. Unscented Kalman filter tutorial. Univ Buffalo Buffalo; 2011. [33] James G, Witten D, Hastie T, Tibshirani R. An introduction to statistical learning: With applications in R. 2013. [34] Cybenko G. Approximation by superpositions of a sigmoidal function. Math Control Signals Systems 1989;2(4):303–14. http://dx.doi.org/10.1007/ BF02551274. [35] ichi Funahashi K. On the approximate realization of continuousmappings by neural networks. Neural Netw 1989;2:183–92. [36] Hagan M, Menhaj M. Training feedforward networks with the marquardt algorithm. IEEE Trans Neural Netw 1994;5(6):989–93. http://dx.doi.org/10.1109/72. 329697. [37] Nocedal J, Wright SJ, editors. Numerical optimization. Springer-Verlag; 1999, http://dx.doi.org/10.1007/b98874. [38] SPS. Sandwich plate systemheavy engineering composite fromintelligent engineering. 2023. [39] Hudson E, Reynolds P. Design and construction of a reconfigurable pedestrian structure. Exp Tech 2016;41. http://dx.doi.org/10.1007/s40799-016-0144-3. [40] AR140. ANSYS user’s manual. 2011. [41] Maia N, Silva J. Theoretical and experimental modal analysis. Engineering dynamics series, Research Studies Press; 1997. [42] Fox RL, Kapoor MP. Rates of change of eigenvalues and eigenvectors. AIAA J 1968;6(12):2426–9. http://dx.doi.org/10.2514/3.5008. [43] Jiménez-Alonso JF, Sáez A, Pavic A, Hudson J. Maximum likelihood methods for finite elementmodel updating of civil engineering structures: A comparative study. In: Proceedings of the 4th international conference on mechanical models in structural engineering cMMoST 2017. 2017, p. 519–32. [44] Kermanshahi B. Design and application of neural networks. Shokodo Tokyo; 1999. [45] Martínez A. Iglesia de san jerónimo (granada), URL http://www.ual.es/ ideimand/iglesia-de-san-jeronimo-granada/. [46] Augenti N, Parisi F, Acconcia E. Mada: online experimental database for mechanical modelling of existing masonry assemblages. 2012. [47] Solutions SV. Artemis modal 5.0. user’s guide. 2015. [48] Wang T, Celik O, Catbas N, Zhang L. A frequency and spatial domain decomposition method for operational strain modal analysis and its application. Eng Struct 2016;27:62–6. http://dx.doi.org/10.1016/j.engstruct.2016.02.011. [49] Peeters B, De Roeck G. Stochastic system identification for operational modal analysis: A review. Trans ASME, J Dyn Syst Meas Control 2001;123. http: //dx.doi.org/10.1115/1.1410370.