A finite element model updating method based on global optimization
Maria Girardi, Cristina Padovani, Daniele Pellegrini, Leonardo Robol
aa r X i v : . [ c s . C E ] J u l A finite element model updating methodbased on global optimization
Maria Girardi a , Cristina Padovani a , Daniele Pellegrini a, ∗ , Leonardo Robol a,b a Institute of Information Science and Technologies “A. Faedo”, ISTI-CNR, Pisa, Italy. b Department of Mathematics, University of Pisa, Pisa, Italy.
Abstract
Finite element model updating of a structure made of linear elastic mate-rials is based on the solution of a minimization problem. The goal is tofind some unknown parameters of the finite element model (elastic moduli,mass densities, constraints and boundary conditions) that minimize an ob-jective function which evaluates the discrepancy between experimental andnumerical dynamic properties. The objective function depends nonlinearlyon the parameters and may have multiple local minimum points. This pa-per presents a numerical method able to find a global minimum point andassess its reliability. The numerical method has been tested on two simu-lated examples – a masonry tower and a domed temple – and validated viaa generic genetic algorithm and a global sensitivity analysis tool. A real casestudy monitored under operational conditions has also been addressed, and ∗ Corresponding author
Email addresses: [email protected] (Maria Girardi), [email protected] (Cristina Padovani), [email protected] (Daniele Pellegrini), [email protected] (Leonardo Robol)
Preprint submitted to Elsevier July 2, 2020 he structure’s experimental modal properties have been used in the modelupdating procedure to estimate the mechanical properties of its constituentmaterials.
Keywords:
Modal analysis, finite elements, model updating, globaloptimization, sensitivity, masonry constructions
1. Introduction
Finite element (FE) model updating is an essential component of numer-ical simulations in structural engineering [24], [29], [41]. It aims to calibratethe FE model of a structure in order to match numerical results with those ob-tained via experimental vibration tests. The calibration allows determiningunknown structure’s characteristics, such as material properties, constraints,and boundary conditions. While the main advantage of such calibration isan updated FE model that can be used to obtain more reliable predictionsregarding the dynamic behaviour of the structure, a further important ap-plication of model updating is damage detection [61], [37], [31].FE model updating consists of solving a constrained minimum problem,the objective function being the distance between experimental and numeri-cal quantities, such as the structure’s natural frequencies and mode shapes.Numerical modal properties depend on some unknown parameters, whichmay suffer from a high degree of uncertainty mainly connected to the lackof information about both the structure’s constituent materials and the in-teractions among its structural elements. In order to reduce the number of2nknown parameters and make the minimum problem more manageable, it ispossible to resort to sensitivity analysis [58], [52], [36], [44], [63], which allowsassessing the influence of the parameters on the modal properties in order toexclude the less influential parameters from the model updating process.Although application of FE model updating to historic masonry build-ings is relatively recent, the literature on the subject is plentiful, [2], [30],[55], [13], [54], [51], [3], [20], [15], [17], [16], [18], [25], [28], [38], [62], [1], [10],[12], [26], [27], and focused on case studies of historical interest for which avibration-based model updating is conducted. Preliminary FE models arecalibrated using the modal properties determined through system identifi-cation techniques. In the majority of the papers cited above the FE modalanalysis is conducted using commercial codes, and the model updating pro-cedure is implemented separately.Many papers have adopted a trial and error approach (see, for example,[20], [13]), in which a manual fine-tuning procedure is used for FE modelupdating. Such an approach is impractical when the number of free parame-ters or the size of the model is large, in which case recourse to an automatedmodel updating becomes more advantageous.The minimum problem stemming from FE model updating, whose ob-jective function may have multiple local minima, can be solved via local orglobal minimisation procedures [57]. The former may be based on trust-region schemes [19], while the latter rely on both deterministic and stochas-tic approaches, which encompass genetic, simulated annealing and particle3warm algorithms.A deterministic approach to the optimisation using multi-start methodsto avoid local minima has been proposed in [27]. In this work the globalminimum point is selected from among several local minima calculated usingdifferent starting points chosen via the Latin Hypercube Sampling (LHS)method [42].A similar approach is adopted in [61] and [8], where the global optimiza-tion technique ”Coupled Local Minimizers”, based on pairwise state synchro-nization constraints, turns out to be more efficient than the multi-start localmethods which rely on independent runs.As far as sensitivity analysis is concerned, several parameter selectionmethods are available for choosing the unknown parameters that should beconsidered in the FE model updating. Most are based on the matrix of lo-cal sensitivities, whose entries usually contain the partial derivatives of thenumerical frequencies calculated at a fixed parameter vector [44]. Local sen-sitivity analysis (LSA) can only provide information about the behaviour ofthe frequencies in a neighbourhood of the given parameter vector and is thusunable to provide any insight into the most relevant parameters influencingthe frequencies. On the other hand, global sensitivity analysis (GSA) [58]provides a global measure of the dependence of the frequencies on the param-eters and represent a preliminary step in the model updating process, whenthe number and influence of the parameters are uncertain. Before tacklingthe optimization problem, it is worth mentioning, by way of example, the4SA applications described in [15] and [27]. In particular, in [15] the resultsof a global sensitivity analysis based on the elementary effect (EE) methodare compared with the results of a local sensitivity analysis, showing that theformer performs better than the latter in model updating of the church of S.Maria del Suffragio in L’Aquila (Italy). Instead, in [27] an average sensitivitymatrix is calculated via the LHS method, which is subsequently adopted tocalibrate the Brivio bridge, a historic concrete structure in Lombardy, Italy.A numerical method for solving the nonlinear least squares problem in-volved in model updating has been proposed in [34] and [35]. The algo-rithm, based on the construction of local parametric reduced-order modelsembedded in a trust-region scheme, was implemented in NOSA-ITACA, anoncommercial FE code developed by the authors [14], [32]. Similar ap-proaches are described in [59] and [27], where the numerical tools expresslydeveloped for model updating are linked to commercial finite element codesused as a black-box within the framework of an iterative process. In par-ticular, [59] presents the MATLAB tool PARIS for automated FE modelupdating. PARIS is a research freeware code linked to the commercial soft-ware SAP2000, which has been applied to full-scale structures for damagedetection purposes. The MATLAB procedure presented in [27] relies insteadon ABAQUS and its efficiency is tested on a historic concrete bridge. Un-like the numerical procedures available in the literature, the algorithm forsolving the constrained minimum problem presented in [34] and [35] takesadvantage of the fact that the NOSA-ITACA source code is at the authors’5isposal. This allows exploiting the structure of the stiffness and mass ma-trices and the fact that only a few of the smallest eigenvalues have to becalculated. To compute these accurately, the natural choice is a (inverse)Lanczos method. When a parametric model is given, the Lanczos projec-tion can be interpreted as a parameter dependent model reduction, wherebyonly the relevant part of the spectrum is matched. The Lanczos projection,combined with a trust-region method, allows matching the experimental fre-quencies with those predicted by the parametric model. This new procedurereduces the overall computation time of the numerical process and turns outto have excellent performance when compared to general-purpose optimizers.In addition, as the procedure described in [34] and [35] allows calculating thesingular value decomposition of the Jacobian of the residual function (thedifference between experimental and numerical dynamic properties) at theminimum point, it makes it possible to assess the reliability of the parameterscalculated and their sensitivity to noisy experimental dynamic properties.In this paper, the numerical method proposed in [34] and [35] to solve theconstrained minimum problem encountered in FE model updating is modifiedin order to calculate a global minimum point of the objective function inthe feasible set. This work is based on a deterministic approach, unlike therelatively recent large body of literature focused on stochastic model updating[40], [63], which aims to take into account and assess the uncertainties in bothexperimental data and numerical models as well.Section 2 recalls the formulation of the optimization problem related to6E model updating. Then the global optimization method integrated intoNOSA-ITACA is described, and some issues related to the reliability of therecovered solution are presented and discussed. In particular, once the op-timal parameter vector has been calculated, we introduce two quantities,which involve the partial derivatives of the numerical frequencies with re-spect to the parameters and provide a measure of how trustworthy the singleparameter is. Section 3 is devoted to testing the numerical method on twosimulated examples: a masonry tower and a domed temple, which highlightthe capabilities and features of the global optimization algorithm proposed inSection 2. For the sake of comparison, we also ran a global optimizer basedon a genetic algorithm available in MATLAB. Such comparisons highlightedthe excellent performance of the proposed method in terms of both compu-tation time and number of evalu ations of the objective function. Section 4presents a real case study, the Matilde donjon in Livorno. This historic tower,which is part of the Fortezza Vecchia (Old Medici Fortress), was subjectedto ambient vibration tests under operational conditions and its experimentaldynamic properties used in the model updating procedure.
2. The numerical method masonry-like materials, which models masonry as an isotropic nonlinear elastic mate-rial with zero or weak tensile strength and infinite or bounded compressivestrength [23], [39]. In recent years, the code has been updated by addingseveral features which now enable it to also to perform modal analysis [4],[5], [6], linear perturbation analysis [50], [33], [49] and model updating [34],[35], [22]. The following subsection 2.1 presents the FE model calibrationas a minimum problem and recalls the algorithm for model updating imple-mented in NOSA–ITACA described in [34] and [35] (to which the reader isreferred for a detailed description). The new features implemented in thecode are explained in detail in subsections 2.2, 2.3 and 2.4.
The term model updating refers to a procedure aimed at calibrating a FEmodel in order to match the experimental and numerical dynamic properties(frequencies and mode shapes) of a structure. It is naturally defined as aninverse problem obtained from modal analysis, which in turn relies on the8olution of the generalized eigenvalue problem Ku = ω Mu , (1)where K and M ∈ R n × n are respectively the stiffness and mass matrices ofthe structure discretized into finite elements, with n the total number of de-grees of freedom. Both K and M are usually sparse and banded, symmetricand positive definite. The eigenvalue ω i is linked to the structure’s frequency f i by the relation f i = ω i / (2 π ), and the eigenvector u ( i ) represents the cor-responding mode shape. The model updating problem can be formulated asan optimization problem by assuming that the stiffness and mass matrices, K and M , are functions of the parameter vector x containing the unknowncharacteristics of the structure (mechanical properties, mass densities, etc.), K = K ( x ) , M = M ( x ) , x ∈ R p . (2)The set Ω of valid choices for the parameters is a p -dimensional boxΩ = [ a , b ] × [ a , b ] ... × [ a p , b p ] , (3)for certain values a i < b i for i =1.... p .By taking (2) into account, equation (1)becomes K ( x ) u ( x ) = ω ( x ) M ( x ) u ( x ) . (4)9he ultimate goal is to determine the optimal value of x that minimizesthe objective function φ ( x ) defined by φ ( x ) = q X i =1 w i [ f i ( x ) − b f i ] (5)within box Ω.The objective function involves the frequencies and therefore dependsnonlinearly on x . We denote by b f the vector of the q experimental frequenciesto match, and by f ( x ) = π p Λ ( x ) the vector of the numerical frequencies,with Λ ( x ) being the one containing the smallest q eigenvalues of Eq. (4),increasingly ordered according to their magnitude. The number p of param-eters to be optimized is expected to be no greater than q . The vector w inEq. (5) encodes the weight that should be given to each frequency in theoptimization scheme. Usually, to obtain relative accuracy on the frequencies w i is chosen equal to b f − i .The objective function φ ( x ) may have several local minima in set Ω,and many numerical methods are available to solve the minimum problem.These include local techniques, often based on trust-region schemes, as well asglobal techniques such as multi-start, genetic and particle swarm algorithms,the last belonging to the class of stochastic methods.A numerical method to find a local minimum point of the objective func-tion introduced above is proposed in [34] and [35], where the authors describea new algorithm based on construction of local parametric reduced-order10odels embedded in a trust-region scheme, along with its implementationinto the FE code NOSA-ITACA. When the FE model depends on parame-ters, as in Eq. (4), and the number n of degrees of freedom is very large, itis convenient to build small-sized, reduced models able to efficiently approxi-mate the behaviour of the original model for all parameter values. Such suchreduced models have been obtained in [34] and [35] through modificationof the Lanczos projection scheme used to compute the first eigenvalues andeigenvectors in Eq. (4) and to create a local model of objective function (5)that is not costly to evaluate and is at least first-order accurate. This localmodel is then used in the region in which it is accurate enough to provide use-ful information on the descent directions; this can be guaranteed by suitablyresizing the trust region, if necessary. It has been be proved that, when thelocal models are accurate, convergence to a local minimizer is guaranteed. Several approaches can be adopted to minimize the objective function(5) in the feasible set Ω. They can be summarized as follows, ordered byincreasing difficulty:1. Find a local minimum point of the objective function in Ω.2. Search for the global minimum point of the objective function in Ω.3. Identify all the local minimum points in Ω and hence, by assuming theyare isolated, recover the global minimum as well.In engineering applications the third approach is the most desirable. Not11nly does it guarantee discovering the most ”likely” parameters, but alsoprovides other values that might be equally acceptable in terms of matchingthe structure’s frequencies. Engineering judgment, something complicatedto insert into an objective function, will then guide the choice of the mostlikely parameter values. In practice, the first approach is easier and alsocomputationally less demanding than both the others, so it is often optedfor.Herein we propose a heuristic strategy to improve the globalization prop-erty of the method introduced in [34] and recalled in the preceding subsec-tion. The goal is to improve the robustness of the method, while partiallyaddressing approaches 2 and 3, without increasing the computational costexcessively. Due to the heuristic nature of the method, from a theoreticalpoint of view, it is impossible to guarantee that all the local minima will befound, but the effectiveness and robustness of the method can be demon-strated through a few practical examples, which are described in the nextsection.The proposed algorithm implemented in NOSA–ITACA code can be sum-marized in the following steps:(a) A local minimum is calculated on the original feasible set Ω = [ a , b ] × . . . × [ a p , b p ], using the method from [34] and assuming the mid-point ofΩ as starting point .(b) For j = 1 , ..., p , let us define m j = ( a j + b j ) and decompose the box Ω12nto the union of 2 p sets of the type¯Ω = I × . . . × I p (6)with I j ∈ { [ a j , m j ] , [ m j , b j ] } , j = 1 , ..., p. (7)(c) A local minimum point is then calculated on each of the subsets definedabove (which have disjoint inner parts), starting at their mid-points. Ifin all the subproblems, the minima coincide with that calculated at step(a), or are on the boundary, then the method stops. Otherwise, therecursion continues on the subsets where new local minima have beenidentified by following the process described in step (b).The method proposed here can run into difficulties when considering alarge number of parameters, as the number of subproblems to solve growsexponentially. However, the following numerical experiments will show thatit is still feasible for several cases of interest.Multi-start optimization approaches are commonly used to find globalminima, for example in [27] the starting points are determined via a LatinHypercube Sampling method and a set of local minimum points found, amongwhich the global minimum point is identified. The algorithm proposed heredoes not execute a fixed number of runs, one for each starting point, but isbased on a recursive procedure, which stops according to a given criterion.Like multi-start methods, the proposed procedure provides a set of local13inimum points, including the global one.The steps laid out above omit one aspect that is rather subtle and requirescareful treatment: how to identify two minimum points. When working infloating-point arithmetic, and using a stopping criterion linked to a specifiedtolerance, two different approximations x and x can be obtained startingfrom two different values for the parameters, even in the case of a singleminimum point. It is therefore essential to be able to distinguish situations inwhich these parameters represent two different minimum points from wheninstead they are just small perturbations of the same minimum point, asexplained in detail in the following subsection. This section is devoted to the open question posed in the foregoing, thatis, how to recognise when two minimum points “coincide”, up to some tol-erance. To answer this question, it is necessary to specify this concept moreclearly. Before addressing this issue, it is worth recalling that the problem ofminimizing function φ in set Ω is a particular inverse problem, as it aims tocalculate the unknown parameters of the FE model of the structure under ex-amination using measurements carried out on it. Analysing minimum pointsprovides a measure of how reliably each parameter has been determined, andcan identify (at the first order) those parameters which only weakly influencethe numerical frequencies, and as such, cannot be reliably determined by theinverse problem. 14ccording to (5) and neglecting vector w for the sake of simplicity, theobjective function under consideration has the form, φ ( x ) = k f ( x ) − b f k , with f ( x ) = f ( x )... f q ( x ) . (8)Let x be a local minimum point of the objective function and assume, upto performing a parameter rescaling, that x is the vector with all componentsequal to 1.Assuming that the objective function is sufficiently regular, the first-orderconditions for x to be a local minimum point imply ∇ φ ( x ) = 0. However,in practical situations vector f is known only approximately, with a tolerance ǫ , so it is possible to introduce a definition of pseudominimum set which isrobust to perturbation.Given x such that ∇ φ ( x ) = 0, we define the ǫ -pseudominimum set at x as follows P ǫ ( φ, x ) = { x | ∃ δ f ∈ R q with k δ f k ≤ ǫ, ∇ φ δ f ( x ) = 0 } , (9)where φ δ f ( x ) = k f ( x ) − b f − δ f k , (10)which is equivalent to considering the set of minimum points of the objectivefunction for close-by frequency configurations, which are acceptable given a15ertain tolerance, ǫ , chosen by the user.In other words, given two local minimum points x and x calculatedvia the scheme described in the foregoing, the two points actually representthe same “numerical” minimum if x ∈ P ǫ ( φ, x ). Note that this relation issymmetric , that is, x ∈ P ǫ ( φ, x ) ⇐⇒ x ∈ P ǫ ( φ, x ), so this definition isconsistent.At first glance, computing P ǫ ( φ, x ) might appear difficult, but by statingthat k x − x k is expected to be small, its determination can be made using afirst-order expansion around x , which makes the problem more manageable.Let us write the first-order expansion of function f ( x ) f ( x ) = f ( x ) + ∇ f ( x )( x − x ) + O ( k x − x k ) , (11)where ∇ f ( x ) denotes the Jacobian of f ( x ) at x = x . In a neighbourhoodof x the Jacobian ∇ f ( x ) can be approximated by the expression ∇ f ( x ) = ∇ f ( x ) + O ( k x − x k ) . (12)From the relation 12 ∇ φ ( x ) = ∇ f ( x ) T ( f ( x ) − b f ) , (13) It is however not transitive, so it does not define an equivalence relation. The dependency of the eigenvalues on the parameters is analytic almost everywherein the domain, hence the Taylor expansion performed here can be rigorously justified. ∇ φ ( x ) = (cid:0) ∇ f ( x ) T + O ( k x − x k ) (cid:1) ( f ( x ) − b f )+ ∇ f ( x ) T ∇ f ( x )( x − x ) + O ( k x − x k ) . (14)We now make the simplifying assumption that the match between theexperimental frequencies and the computed ones is good enough for the termin the first line of (14) to be negligible, and hence12 ∇ φ ( x ) . = ∇ f ( x ) T ∇ f ( x )( x − x ) . (15)Then, from (14), taking into account that k δ f k ≤ ǫ , we get an analogueexpression for ∇ φ δ f ( x ),12 ∇ φ δ f ( x ) . = ∇ f ( x ) T ∇ f ( x )( x − x )) − ∇ f ( x ) T δ f . (16)For the sake of simplicity, let us denote with the same symbol, P ǫ ( φ, x ),the set computed by replacing ∇ φ δ f ( x ) with the above approximation. Then,it follows that P ǫ ( φ, x ) = (cid:8) x | ∃k δ f k ≤ ǫ, ∇ f ( x ) T ∇ f ( x )( x − x ) = ∇ f ( x ) T δ f (cid:9) . (17)Let UΣV T = ∇ f ( x ) T be the singular value decomposition (SVD) of17 f ( x ) T . By virtue of the fact that δ f is arbitrary, and the multiplicationby unitary matrices leaves the Euclidean norm unchanged, it is possible torewrite the set in (17) as follows P ǫ ( φ, x ) = (cid:8) x | k ΣU T ( x − x ) k ≤ ǫ (cid:9) . (18)A SVD can be compute with O ( q p ) flops, assuming q ≥ p , and is thereforea negligible cost in the proposed algorithm. Note in particular that the costof computing this set is independent of n , the degrees of freedom in the FEmodel. Hence, (18) is easily verifiable in practice, and has been implementedas a test in the algorithm described in the foregoing. The algorithm returnsthe matrices Σ and U , which can be used to construct the ellipsoid P ǫ ( φ, x ),which describes, at the first-order, the level of accuracy attained in the spaceof parameters. In addition, the SVD of the Jacobian can be used to compute,for each parameter x j , the quantities ζ j and η j , as described in the nextsubsection. Generally, experimental frequencies may not be accurate, since they arederived by analyzing measured data that may be contaminated by environ-mental noise. Thus, when minimizing objective function (5), one has toensure that the optimal parameters are well-defined and robust to perturba-tions in the data b f .This analysis is only relevant in a neighbourhood of the minimum point:18he behaviour of the objective function elsewhere does not influence the con-ditioning of the optimization problem.A complete description of the parameters space and the directions wherethe problem is well- or ill-defined can be given by computing the SVD of theJacobian, as is widely referenced in the numerical optimization literature andpointed out for the problem at hand in [35]. Nevertheless, if the dimension ofthe parameter space is greater than three, giving a meaningful interpretationto these directions can be difficult; hence, we introduce two quantities whichare easier to interpret and convey the same information.Let b x be a local minimum point of the nonlinear objective function (5).We assume that function f ( x ) has been properly scaled so that both b x and b f are vectors of all ones, and we replace f ( x ) with its first-order expansion(11), having set x = b x . We may now define the following parameters foreach j = 1 , . . . , p ζ j := (cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂x j (cid:13)(cid:13)(cid:13)(cid:13) . η j := min v ∈S j (cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ v (cid:13)(cid:13)(cid:13)(cid:13) , (19)where ∂ f ∂ v denotes the directional derivative, and set S j is defined as follows S j := v v ∈ R p | (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) v v (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ , v ∈ R j − , v ∈ R p − j . (20)Note that set S j contains, in particular, the j -th vector e j of the canonical19asis of R p , and therefore it must hold that η j ≤ ζ j . Intuitively, S j is theset of directions where j − th parameter is forced to change at “unit speed”,while the others can change at some other speed, but are still bounded inthe Euclidean norm by 1. Taking the minimum of the directional derivativesin S j is equivalent to finding the direction in the parameter space with theslowest growth of f ( x ), in which parameter x j is involved.Hence, we can make the following remarks: • If η j is small (i.e., η j ≪ x j is forced to change, but f ( x ) varies slowly; hence, determination of x j might be subject to noise. If, on the other hand, η j ≫
0, then itsdetermination through the optimization problem is robust to noise. • If ζ j is small, then when x j changes, the frequencies are nearly unaf-fected; hence, there is no information on x j that can be obtained bysolving the optimization problem. On the other hand, if ζ j is large,then it cannot be guaranteed that x j is not affected by noise, but thereis at least one direction in the parameter space involving x j that canbe reliably determined.The direction mentioned above can be determined from the SVD of theJacobian ∇ f ( b x ) = UΣV T , as described in [35]. However, parameters ζ j and η j are easier to read, and we have the following trichotomy:(i) η j ≤ ζ j ≪
1: parameter x j cannot be reliably determined, as no infor-mation on it is encoded in the optimization problem.20ii) 0 ≪ η j ≤ ζ j : parameter x j can be reliably determined from the data,even if it is subject to noise. The amount of noise that can be toleratedis bounded in norm by η j .(iii) η j ≪
1, but ζ j ≫
0: there is some information on parameter x j en-coded in the problem, but the result will not be free of noise. To findthe directions which can be “trusted”, one has to look at the right sin-gular vectors corresponding to large singular values in the SVD of theJacobian.It is immediately clear that ζ j can be computed directly by taking thenorms of the columns of the Jacobian. Computing η j , on the other hand,requires some more effort. Let us temporarily drop the requirement that k [ v T v T ] k < v can be found by solving anunconstrained linear least square problem, and in particular we have v = v v , with v v = −∇ f ( b x ) † j ∇ f ( b x ) e j , (21)where ∇ f ( b x ) j is the Jacobian without the j -th column, and the symbol † denotes the Moore-Penrose pseudoinverse. If k [ v T v T ] k is less than 1,then v in (21) is the minimizer for the constrained problem in (19) as well.Otherwise, an explicit formula is not available and we use the orthogonalprojection of the computed v onto S j as a starting point and determine the21olution by solving a constrained nonlinear least square problem. For solutionof this problem, we rely on the SQP algorithm described in Chapter 18 of[46].
3. Application to simulated case studies
In order to test the method described in section 2, two artificial exampleshave been proposed. In both cases, the structure’s free parameters are as-signed, and a preliminary numerical modal analysis is performed to evaluatethe corresponding frequencies and mode shapes. Subsequently, the numer-ical frequencies are employed as input to the model updating procedure torecover the original parameters. The first example highlights the ability ofthe NOSA–ITACA code to discover more minimum points as compared to ageneric genetic algorithm used to solve the same problem, which is unable tofind more than one point. The second example shows some of the code’s fea-tures, which can help users to choose the most suitable optimal parameterscharacterized by the greatest reliability.The tests, conducted with NOSA-ITACA and MATLAB R2018b, wererun on a computer with an Intel Core i7-8700 running at 3.20 GHz, with64GB of RAM clocked at 2133MHz.The weight vector w is always chosen to be w i = b f − i , which ensuresrelative accuracy of the recovered frequency.22 .1. A masonry tower As a first example, we considered the tower shown in Figure 1. The20 m-high structure has a rectangular cross section of 5 m ×
10 m and wallsof 1 m constant thickness. The tower, clamped at its base, is discretizedinto 2080 eight–node quadrilateral thin shell elements (element number 5 ofthe NOSA-ITACA library [14]) for a total of 6344 nodes and 25376 degreesof freedom. A preliminary modal analysis is performed to evaluate the fre-quencies and mode shapes under the assumptions that the tower is made of ahomogeneous material with Young’s moduli E = E = 3 .
00 GPa (see Figure1), Poisson’s ratio ν = 0 . ρ = 1835 . / m . The vectorof the corresponding natural frequencies obtained with the above parametersis b f = [2.670, 4.737, 6.571] Hz . (22)Figure 1 shows the mode shapes corresponding to the first three tower’sfrequencies: the first two modes are bending movements along X and Yrespectively, while the third is a torsional mode shape.23 igure 1: The masonry tower: geometry (length in meters); model created by NOSA-ITACA code; the first three mode shapes. The algorithm described in this paper is used to determine the Young’smoduli E and E of the structure. Putting x = [ E , E ], with the parametervarying within the interval1 .
00 GPa ≤ E , E ≤ .
00 GPa , (23)model updating is conducted considering frequencies b f and b f in case (a),and b f , b f and b f in case (b).The same problems are also addressed with a generic genetic algorithm24denoted by GA) available in MATLAB R2017b, using NOSA–ITACA asa black box, with the aim of comparing the results of the two approachesand test the reliability and robustness of the numerical procedure proposed.Table 1 summarizes the results related to case (a). Note firstly that NOSA–ITACA code finds two minimum points, which correspond to the exact valuesof the known frequencies, while the genetic algorithm calculates only oneminimum, which is expected be the global minimum point. The existence oftwo minimum points is shown in Figure 2, where the plot of the objectivefunction φ ( x ) defined in Eq. (5) is reported in log–scale, as the two elasticmoduli vary. Regarding computation times and the number of evaluationsof the objective function, the numerical procedure implemented in NOSA–ITACA appears to be much more efficient.NOSA–ITACA GAMinimum 1 [3.00; 3.00] GPa [3.02; 2.95] GPaFrequencies [2.670, 4.737] Hz [2.671, 4.732] HzMinimum 2 [4.49; 1.34] GPa –Frequencies [2.670, 4.737] Hz –Computation time 11.50 s 465.03 sNumber of evaluations 41 2600 Table 1: Case (a) – Optimization results, two frequencies and two parameters. igure 2: Case (a) – Objective function vs. E and E . Regarding case (b), the results summarized in table 2 clearly show thesuperior performance of the NOSA–ITACA code in terms of both computa-tion time and accuracy. Figure 3 shows the plot of the objective function φ ( x ), which in this case exhibits one global minimum point.NOSA–ITACA GAMinimum 1 [3.00; 3.00] GPa [3.00; 2.99] GPaFrequencies [2.670, 4.737, 6.571] Hz [2.670, 4.737, 6.571] HzComputation time 7.72 s 497.63 sNumber of evaluations 27 2600 Table 2: Case (b) – Optimization results, three frequencies and two parameters. igure 3: Case (b) – Objective function vs. E and E . Table 3 shows, for each minimum point of cases (a) and (b), the param-eters values ζ j and η j defined in subsection 2.4. In all cases, 0 ≪ η j ≪ ζ j ,which means that every parameter E j has been determined reliably (as isevident in tables 1 and 2) from the data, even if subject to noise. The tablealso report ζ − j and η − j , quantities which provide an estimate of the orderof magnitude of the minimum and maximum percentage error (at the first-order) inherent in estimating the parameters under the hypothesis of a 1%error in the assessment of the experimental frequencies. From the table it isclear that, in the worst-case scenario, parameter estimation will be affected,27t most, by a 6.2% error in both cases (a) and (b).Case Minimum x j ζ j η j ζ − j η − j (a) 1 E E E E E E Table 3: Parameters ζ j and η j for the cases (a) and (b). Let us now consider the domed temple, depicted in Figure 4, consistingof a 5 m high octagonal shaped cloister vault resting on a drum inscribedon a 10 m ×
11 m rectangle. The structure, clamped at its base, is made of4 different materials (Figure 5): material 1 for the dome (orange), material2 for the upper part of the drum (cyan), material 3 for the bottom part ofthe drum (violet) and material 4 for the columns (green). The finite elementmodel, shown in Figure 5, is composed of 31052 hexahedron brick elementsand 41245 nodes for a total number of 123735 degrees of freedom.28 igure 4: Geometry of the domed temple (length in meters). igure 5: Domed temple, mesh and materials. Each color corresponds to a differentmaterial, orange (1), cyan (2), violet (3) and green (4). A preliminary modal analysis is performed to evaluate the structure’sfrequencies assuming the material properties reported in table 4. The vectorof the first eight natural frequencies is b f = [2.19, 2.23, 3.76, 3.83, 4.32, 4.60, 4.72, 8.26] Hz . (24)30aterial Temple portion ρ [kg/m ] E [GPa] ν Table 4: Values of the material properties.
The optimization code implemented in NOSA–ITACA and a generic ge-netic algorithm were run setting x = [ E , ρ , E , E , ρ , E , ρ ], with the fol-lowing bounds 2 .
00 GPa ≤ E j ≤ .
00 GPa , j = 1 , ..., , (25)1600 . / m ≤ ρ j ≤ . / m , j = 1 , , . (26)This choice leaves seven parameters to be optimized, with the sole exceptionof ρ , which was set to the fixed value reported in table 4. Tables 5 and6 summarize the results obtained by NOSA–ITACA code and the geneticalgorithm in terms of optimal parameter values, frequencies, relative errors | ∆ x j | and | ∆ f | , computation time and number of evaluations of the objectivefunction. 31eal value NOSA–ITACA | ∆ x j | [%] GA | ∆ x j | [%] E [GPa] 3.000 2.996 0.13 4.1431 38.10 ρ [kg/m ] 1800.0 1908.9 6.05 1988.6 10.47 E [GPa] 3.500 4.085 16.72 4.0335 15.24 E [GPa] 4.000 4.177 4.43 3.8357 4.11 ρ [kg/m ] 2000.0 2115.9 5.80 2340.1 17.00 E [GPa] 5.000 5.132 2.63 5.6213 12.43 ρ [kg/m ] 2200.0 2272.7 3.30 2397.8 9.00Computation time [s] 14019 103250Number of evaluations 671 10500 Table 5: Optimal parameter values calculated by NOSA–ITACA code and a genetic algo-rithm.
Real value NOSA–ITACA | ∆ f | [%] GA | ∆ f | [%] f [Hz] 2.19 2.18 0.46 2.18 0.46 f [Hz] 2.23 2.22 0.45 2.22 0.45 f [Hz] 3.76 3.75 0.27 3.77 0.27 f [Hz] 3.83 3.83 0.00 3.83 0.00 f [Hz] 4.32 4.31 0.23 4.31 0.23 f [Hz] 4.60 4.60 0.00 4.61 0.22 f [Hz] 4.72 4.72 0.00 4.72 0.00 f [Hz] 8.26 8.25 0.12 8.24 0.24 Table 6: Frequencies values corresponding to the parameters’ optimal values recovered byNOSA–ITACA code and a genetic algorithm.
The results above highlight that: (i) the numerical procedure imple-mented in NOSA–ITACA is less time–consuming than the genetic algorithm,the computation time of the former being ten times lower than that of thelatter; (ii) the optimal values of the Young’s moduli calculated by NOSA–ITACA are affected by a maximum relative error of 17%, against 38% ofthe genetic algorithm; (iii) the maximum relative error on mass density is32bout 6% for NOSA–ITACA and 17% for the genetic algorithm; (iv) eventhough the optimal value of some mechanical characteristics is affected byhigh error, the maximum relative error on the frequencies is about 0 .
5% forboth numerical methods.To investigate the robustness and reliability of the solution found, theparameters values ζ j and η j defined in subsection 2.4 are reported in table 7with their respective inverse values and the relative error | ∆ x j | calculated intable 5. ζ j η j ζ − j η − j | ∆ x j | [%] E · − · − ρ · − · − E · − · − E · − · − ρ · − · − E · − ρ · − Table 7: Parameters ζ j and η j calculated by NOSA–ITACA. The above table shows that the Young’s moduli of materials 1 and 2 (thedome and the upper part of the drum) seem to be irrelevant in the opti-mization process. This fact can be explained by observing the mode shapesrelated to the first eight frequencies, which mainly involve displacement ofthe pillars. It is also interesting to note that the objective function is moreheavily influenced by the dome’s mass density than by its elastic modulus( ζ = 5 . · − versus ζ = 1 . · − ), in line with the fact that thedynamic behavior of the structure is comparable to a cantilever beam with a33ass concentrated at the free end. The Young’s moduli and mass density ofmaterials 3 and 4 seem more reliable than the others, as shown by the valuesof ζ j and η j . Finally, note that the relative error | ∆ x j | made in estimatingthe optimal values of the parameters is always close to the range defined by ζ − j and η − j (at the first-order, under the hypothesis of a maximum error of1% in the assessment of the experimental frequencies).Further information can be achieved by calculating, at the minimumpoint, the scaled Jacobian matrix described in subsection 2.4, . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − . · − − . · − . · − . · − − . · − . · − − . · − (27)The numbers reported in the first three columns of the matrix confirmsthat the temple’s frequencies are weakly dependent on materials 1 and 2.Restricting the attention to the last two columns in matrix (27) (containingthe partial derivatives of the frequencies with respect to E and ρ ) furnishesmore information about the minimum point. The SVD of the restricted34atrix yields the results summarized in table 8, with the singular values σ > σ reported in the first columns, and the corresponding right singularvectors in the second and third columns. The objective function is expectedto have a direction with a weaker influence on the frequencies parallel to z (2) (with constant ratio E / ρ ), which corresponds to the smallest singular value σ = 2 . · − . σ z (1) z (2) · − -6.8971 · − · − · − -7.2408 · − Table 8: Singular values and right singular vectors of the scaled restricted Jacobian matrix.
To investigate how variation in the input (Young’s moduli and the massdensities of the domed temple’s four constituent materials) influence the out-put of the numerical model (the natural frequencies), and thereby test thesensitivity analysis implemented in the NOSA–ITACA code, a Global Sensi-tivity Analysis (GSA) has been performed through the SAFE Toolbox [52],[45] and [53].The SAFE Toolbox, an open–source code implemented in MATLAB, canbe easily linked to simulation models running outside the MATLAB environ-ment, such as the NOSA-ITACA code in the example at hand. The Elemen-tary Effects Test (EET method [43]) is used to evaluate the sensitivity indicesassuming that the eight input parameters (Young’s moduli and the mass den-sities of the four materials) have a uniform probability distribution function,and adopting the Latin Hypercube method [42] as sampling strategy. From35igure 6, where the sensitivity indices calculated via the EET method areplotted, it is possible to deduce that the Young’s moduli of materials 3 and 4affect the numerical frequencies much more than the remaining parameters.These results confirm the information recovered by the quantities ζ j and η j calculated by NOSA–ITACA and reported in table 7. Figure 6: EET sensitivity indices for the first nine frequencies and eight parameters.
Figure 7 shows the two-dimensional scatter plots related to the first nat-ural frequency, which allows identifying pairwise interactions among inputfactors. From the plot, it is clear that no interactions occur between any ofthe parameters, except for the mass density and Young’s modulus of material36: the first frequency depends on the ratio E / ρ (as indicated by the dashedblack arrow).Moreover, the scatter plots in the sixth column highlight the fact thatfor fixed E , the frequency does not depend on other parameters. Similarconsiderations can be made by constructing two-dimensional scatter plots forthe other frequencies. All these considerations corroborate the conclusionsreached through the sensitivity analysis implemented in the NOSA–ITACAcode. It is also worth noting that the computational cost of such a globalsensitivity analysis is very high (Figures 6 and 7 are the results of 1260 FEmodal analysis runs) with respect to the cost of the minimization procedureimplemented in NOSA-ITACA, which provides both the global minimumpoint and an assessment of its reliability.37 igure 7: Two-dimensional scatter plot related to first frequency.
4. Application to a real example: the Matilde donjon in Livorno
The Matilde donjon is a fortified keep belonging to the Fortezza Vecchia(Old Fortress), near the ancient Medici Port of Livorno, Italy (Figure 8).38
The 26 m–high cylindrical tower shown in Figures 9 and 10 has a cross-section with a mean outer radius of 6 m andwalls of 2 . Figure 9: The Matilde donjon (view 1, 2). igure 10: The Matilde donjon (view 3). In October 2017, an ambient vibration monitoring experiment was car-ried out on the tower (see Figure 11, 12, 13). The ambient vibrations weremonitored for a few hours via SARA SS20 seismometric stations (made avail-able by INGV of Arezzo) arranged in different layouts. During the five tests(T1 to T5), each lasting about thirty minutes, two sensors were kept in afixed position–, one at the base (level -2) and the other on the roof terrace(level 2)–, while the remaining sensors were moved to different positions alongthe tower’s height and surrounding area in order to obtain information onthe mode shapes and degree of connection between the Old Fortress’ struc-tures and the tower itself. The sampling rate was set at 100 Hz. All datarecorded have been divided into short sequences, each lasting 1000 seconds(a time window greater than the structure’s fundamental period estimatedby preliminary FE modal analysis), and processed by two different opera-41ional modal analysis (OMA) techniques, through which the tower’ modalparameters were estimated: the Stochastic Subspace Identification covari-ance driven method (SSI–cov) [47] implemented in MACEC code [56] andthe Enhanced Frequency Domain Decomposition method (EFDD) [11] im-plemented by ISTI in Trudi code [48].
Figure 11: Transverse sections of the tower. igure 12: Sensor layout October 2017 – test T1, T2, T3. igure 13: Sensor layout October 2017 – test T4, T5. In total, six vibration modes were identified in the frequency range of2-13 Hz. Table 9 summarizes the results in terms of natural frequencies f ,damping ratios ξ , and MAC values [60] calculated between the correspondingmode shapes estimated via the two OMA techniques.For the sake of brevity, the values shown in the tables correspond to theaverage values of the estimated parameters during each test, all of which arecharacterized by a MPC value [60] greater than 0 . SSI-cov [Hz] ξ SSI-cov [%] f EFDD [Hz] ξ EFDD [%] MAC
SSI-ref,EFDD
Mode 1 2.68 3.47 2.69 2.97 0.99Mode 2 3.37 3.90 3.35 4.11 0.99Mode 3 6.21 1.44 – – –Mode 4 8.10 4.63 8.15 1.14 0.97Mode 5 10.04 5.69 10.06 – 0.97Mode 6 11.95 1.15 12.24 – 0.99
Table 9: Modal parameters of the tower, October 2017.
The two first mode shapes are bending mode along the west-east directionand north–south direction, respectively, while the third mode corresponds totorsional movement of the tower and a deflection of the two lateral wallsconnected to its south–west portion. The other experimental mode shapesare more uncertain: the fourth one is likely a torsion mode shape mixed withbending along north-east/south-west direction, and the fifth and sixth arehigher–order bending mode shapes.
In this subsection, the procedure described in Section 2 is applied to theMatilde donjon. The FE mesh of the tower, shown in Figure 14, consists of52560 isoparametric eight-node brick elements and 64380 nodes, for a totalof 193140 degrees of freedom. The model, as shown in the Figure, includesa portion of the surrounding walls. The bases of the tower and lateral wallsare fixed, and the ends of the walls are prevented from moving along the Xand Y directions. 45 igure 14: FE model of the Matilde donjon.
The numerical procedure has been used to estimate the values of theYoung’s modulus of the inner and outer layers ( E t,i = E t,e = E t ) of thetower’s walls, and Young’s moduli ( E m,i ) of the masonry constituting theFortress’ walls (Figure 15), with x = [ E t , E m, , E m, , E m, ]. These parame-ters have been allowed to vary within the intervals1 .
00 GPa ≤ E t ≤ .
00 GPa , (28)46 .
00 GPa ≤ E m, , E m, , E m, ≤ .
00 GPa . (29) Figure 15: Designated tower materials.
The Poisson’s ratio of masonry is fixed at 0 .
2, the mass density of thetower’s walls is fixed at 1800 kg / m and 2000 kg / m for the inner and outerlayer, respectively, and the mass density of the side walls is taken to be2000 kg / m . The experimental frequencies estimated by the SSI–cov methodare used in the optimization process, hence47 f = [2.68, 3.37, 6.21, 8.10, 10.04, 11.95] Hz . (30)The optimal parameters are reported in table 10: the values of ζ and η guarantee the reliability of E t and E m, , while the constituent materials theremaining walls are marked by uncertainty. The total computation time forthe model updating procedure was 8468.9 s, and the number of evaluations131. x j ζ j η j ζ − j η − j E t [GPa] 2.152 1.627 1.557 0.615 0.642 E m, [GPa] 5.808 9.577 · − · − E m, [GPa] 5.532 6.409 · − · − E m, [GPa] 2.095 6.845 · − · − Table 10: Optimal parameter values calculated by NOSA–ITACA.
Table 11 summarizes the numerical frequencies of the tower correspondingto the optimal parameters and their relative errors | ∆ f | with respect to theexperimental counterparts; | ∆ f | varies between 2 and 3%, except for thethird and sixth frequencies. 48 f i [Hz] f i [Hz] | ∆ f | [%]mode 1 2 .
68 2 .
76 2 . .
37 3 .
33 1 . .
21 6 .
51 4 . .
10 7 .
90 2 . .
04 9 .
81 2 . .
95 11 .
10 7 . Table 11: Experimental frequencies b f and numerical frequencies f calculated for the opti-mal values of the parameters recovered by NOSA–ITACA. As for the simulated example, a GSA has been performed to validatethe results of the sensitivity analysis achieved by NOSA–ITACA. The EETmethod is used to evaluate the sensitivity indices assuming a uniform proba-bility distribution function, for the nine input factors (Young’s modulus andmass density of each material), and the Latin Hypercube as sampling strat-egy; 500 FE modal analyses were carried out. Figure 16 shows that the elasticmoduli of the tower and wall 1 strongly influence the frequency variation ascompared to the others. In particular, the tower’s Young’s modulus impactsall frequencies except for the third, which is instead heavily affected by elasticmodulus E m, , as confirmed by the experimental mode shape which exhibitsa large displacement component corresponding to an out-of-plane deflectionof the wall. The GSA analysis confirms the reliability of the NOSA–ITACAresults. 49 igure 16: EET sensitivity indices for the first sixth frequencies and nine parameters.
5. Conclusions
The present paper proposes an improved numerical method to solve theconstrained minimum problem encountered in FE model updating and cal-culate a global minimum point of the objective function in the feasible set.The global optimization method, consisting of a recursive procedure basedon construction of local parametric reduced-order models embedded in atrust-region scheme, is integrated into the FE code NOSA-ITACA, a soft-ware developed in house by the authors. Along with the global optimizationmethod, some issues related to the reliability of the recovered solution are50resented and discussed. In particular, once the optimal parameter vectorhas been calculated, two quantities involving the Jacobian of the numericalfrequencies provide a measure of how trustworthy the single parameter is.The numerical method has been tested on two simulated examples a masonrytower and a domed temple in order to highlight the capabilities and featuresof the proposed global optimization algorithm. The results of the test cases,validated via a generic genetic algorithm and a global sensitivity analysis,prove the method’s efficiency and robustness. The objective function mayhave multiple local minimum points, and the first example highlights thatthe proposed procedure, unlike a genetic algorithm, can provide a set of localminimum points, including the global one. The second example shows somefeatures of the code, which can help users to choose the most suitable optimalparameters characterized by higher reliability. Comparison of the computa-tion time and number of objective function evaluations highlights that theNOSA-ITACA code performs better than the genetic algorithm. Regardinghow the parameter variations can influence the frequencies of the FE model,the numerical method seems to provide the same information given by aglobal sensitivity analysis. Finally, the paper has addressed a real case studythe Matilde donjon in Livorno. The experimental dynamic properties of thehistoric tower monitored under operational conditions were used in the modelupdating procedure to estimate the mechanical properties of its constituentmaterials. 51 cknowledgements
The authors wish to thank Dr. Riccardo Mario Azzara for having madeavailable the seismic instrumentation used in the experimental tests per-formed on the Matilde donjon.