PhysEx: A novel procedure to estimate full-rank physical matrices of a structure from an incomplete modal model

Journal of Sound and Vibration 510 (2021) 116277 Available online 11 June 2021 0022-460X/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi PhysEx: A novel procedure to estimate full-rank physical matrices of a structure from an incomplete modal model Alvaro Magdaleno∗, Antolin Lorenzana ITAP. Universidad de Valladolid, Paseo del Cauce 59, 47011 Valladolid, Spain ARTICLE INFO Keywords: Dynamic modelling Incomplete model Physical model Structural parameter extraction ABSTRACT This work is devoted to presenting and applying a novel methodology that estimates a linear physical model, represented by full-rank mass, stiffness and damping matrices, capable of replicating the dynamic behaviour of an incomplete modal model. No hypothesis is made with respect to the damping matrix, so the procedure is applicable assuming any linear viscous damping model. It is shown that, if the number of considered modes satisfies a certain restriction with respect to the number of degrees of freedom, the problem has infinite solutions that differ in their exogenous eigenvalues. By controlling their effect on the dynamic behaviour of the model in the frequency range of interest, a proper final solution is obtained. The methodology is applied to a three-degrees-of-freedom discrete model and a discretized cantilever beam. Comparisons are made in both frequency and time domains which reveal that the methodology provides physical models that accurately replicate the dynamic behaviour of the incomplete modal model. Even for complete models, where the solvability condition is not met, the methodology manages to provide meaningful results in some situations. To prove the usefulness of the estimated models, a mass is added on one DOF of both examples to show that they reproduce the same effects as the original modified ones. Finally, in the interests of reproducibility, all models and algorithms presented in this article have been made publicly available. 1. Introduction One desirable requirement in structural engineering is the ability to build up computational models capable of accounting for the behaviour of a certain structure as close as possible to reality. This enables the engineer to refine its design, accurately predicting its response when different kinds of loads are applied, specially those that cannot be tested, such as earthquakes. In addition, the structural behaviour can be simulated in different scenarios, such as the structure undergoing some kind of structural modification, damage or having a vibration mitigation device incorporated. When the dynamic behaviour is of interest, it is hard to accurately predict the structural properties before building. Only a strictly balanced amount of mass, stiffness and damping of the structural elements can provide a realistic dynamic response and, even if a highly detailed model is elaborated taking into account all the structural and non-structural sources of mass and stiffness, thus leading to huge models that are very cumbersome to work with, it is still very hard to accurately determine the damping contribution of each structural element. Therefore, it is a common practice for the final computational model, truly representing the dynamic behaviour of the structure, to be issued after having carried out a series of dynamical tests to the real, already built structure. ∗Corresponding author. E-mail addresses: [email protected] (A. Magdaleno), [email protected] (A. Lorenzana). https://doi.org/10.1016/j.jsv.2021.116277 Received 13 June 2020; Received in revised form 26 April 2021; Accepted 8 June 2021 Journal of Sound and Vibration 510 (2021) 116277 2 A. Magdaleno and A. Lorenzana The relationship between the experimental data and the computation model divides the main existing techniques into two large groups: those which aim to improve or update an already existing computational model to fit the data; and those intended to estimate a new model from the data. The first group of algorithms, usually referred to as Model Updating techniques [1–3], normally provides good results after a potentially long iterative updating process. The final results depend on the quality and size of the starting model, which must contain a sufficient number of adjustable parameters to fully adapt its response to the experimental data [4,5]. On the other hand, the techniques whose purpose is to directly estimate a model from the experimental data, commonly known as Model Parameter Extraction techniques, aim to develop a set of mathematical entities that, combined, reproduce the behaviour of the structure in the frequency range of interest. One example of this, assuming a linear behaviour, is the modal model, which, by means of a set of natural frequencies, damping ratios and scaled mode shapes, is able to predict the time response of the structure for any kind of applied load [6,7]. The goodness of the estimation, in this case, depends on the quality of the measurements and the ability to extract the modal model from the experimental data. In addition, although suitable for simulating a wide range of actions, some typical engineering tasks, such as structural modification or vibration mitigation devices tuning, may be difficult to undertake through purely mechanical approaches. This work aims to go one step further to estimate a linear physical model from a given modal model. A physical model is composed of mass, stiffness and damping matrices, similar to those that can be assembled through any other computational approach, such as the Finite Element Method. Being issued from a modal model and not from the assembly of smaller matrices associated to specific structural elements, they can potentially have any internal form, thus giving them the ability to better adapt the measured dynamics. The problem can first be addressed by considering the same number of modes 𝑚as measured degrees of freedom 𝑛 (complete models), as in [8], but this is not often the case. Assuming that a sufficient number of sensors have been used to measure the modal properties of the structure, spatial incompleteness cannot occur [9]. It is common for the number of identified modes to be fewer than the number of monitored degrees of freedom, and it is necessary to deal with incomplete modal models. Berman [10,11] warned about the solution to this problem not being unique and some authors [12,13] have explored it since then by inverting, or pseudo-inverting, the orthogonality properties of the mode shapes. This approach can easily be applied if viscous proportional damping is assumed, since it eases the formulation thanks to the real-valued mode shapes. However, some authors addressing complex-valued mode shapes proposed transforming them into real-valued ones [14] or choosing another kind of approach [15,16]. In fact, the damping matrix has traditionally been the most complicated one to estimate and most recent works in this topic have focused on this [17–21]. In these works, both the mass and stiffness matrices are assumed to be previously known or estimated by other means, such as a finite element model. In this paper, the problem of estimating the three matrices of a physical model from an incomplete complex-valued modal model is addressed, so a general damping model is assumed. The methodology, PhysEx (named after the terms PHYSical parameter EXtraction), is developed by analysing the different requirements that the physical matrices must meet in order to reproduce the objective dynamic behaviour. These requirements are coded as a set of linear equations that must be solved for the elements of the physical matrices. It is demonstrated that, if the number of considered modes is fewer than a fraction of degrees of freedom, the problem does not have a unique solution. In this case, the infinite possible solutions differ in the so-called exogenous modes, i.e., the modes that mathematically belong to the physical model which are not the objective modes to reproduce. By controlling their effect on the dynamic behaviour in the frequency range of interest, a final unique solution can be calculated through an optimization process. As a result, a set of full rank physical matrices is obtained, even if the initial modal model is incomplete. As will be seen in the corresponding section, although the terms which compose the matrices do not have any direct physical meaning, the estimated physical models are useful to perform some tasks. As a proof, a certain amount of mass is added to the estimated mass matrix to show that this reproduces the same effects as the addition of the same amount of mass to the original model. Future works will be devoted to covering different structural modifications, such as stiffness and damping modifications, as well as the addition and tuning of inertial mitigation devices. The article is structured as follows: Section 2presents the essential modal analysis background. The methodology is developed in Sections 3and 4gives two application examples consisting in a three-degrees-of-freedom discrete model and a discretized cantilever beam. A dataset containing all the models and code developed in MATLAB for this work is available in [22] for reproducibility purposes. Finally, Section 5ends the article by summing up the main conclusions. 2. Theoretical background The linear behaviour of a discretized structure can be described using its mass (𝐌), damping (𝐂) and stiffness (𝐊) matrices as shown in Eq. (1), where 𝐪(𝑡)and 𝐅(𝑡)are column vectors standing for, respectively, the 𝑛degrees of freedom (DOF) used to describe the motion of the structure and the external forces acting on them. The dot operand () stands for a time derivative. 𝐌 𝐪(𝑡) + 𝐂 𝐪(𝑡) + 𝐊𝐪(𝑡) = 𝐅(𝑡)(1) When a non-proportional linear damping model is assumed, the system of equations in Eq. (1) cannot be directly uncoupled through its modal properties [6]. One way to proceed is by transforming the equations to the state space, by means of the state vector 𝐱(𝑡)=[𝐪𝑡(𝑡) 𝐪𝑡(𝑡)]𝑡and an input vector 𝐮(𝑡)=[𝐅𝑡(𝑡) ∅]𝑡, leading to Eq. (2), where ∅stands for a matrix of zeros of a proper dimension. 𝐀 𝐱(𝑡) + 𝐁𝐱(𝑡) = 𝐮(𝑡) 𝐀=[𝐂 𝐌 𝐌∅]𝐁=[𝐊∅ ∅ −𝐌](2) Journal of Sound and Vibration 510 (2021) 116277 3 A. Magdaleno and A. Lorenzana The system of equations in Eq. (2) can be uncoupled by means of its 𝑛complex conjugate pairs of eigenvectors ( 𝜽𝑟, 𝜽∗ 𝑟), where 𝑟= 1..𝑛 and ∗stands for the complex conjugation. They can be calculated, together with their corresponding eigenvalues (𝑠𝑟, 𝑠∗ 𝑟), by solving the eigenproblem in Eq. (3). For the viscous damping model, the eigenvalues have the form 𝑠𝑟= −𝜔𝑟𝜁𝑟+ j𝜔𝑟√1 − 𝜁2 𝑟, where jstands for the imaginary unit, 𝜔𝑟is the natural frequency (in rad/s) and 𝜁𝑟is the viscous damping ratio. (𝐀𝑠𝑟+𝐁) 𝜽𝑟= 0 (3) It is worth noting that, due to the form of the state vector 𝐱(𝑡), the eigenvectors  𝜽𝑟are composed of the mode shapes (𝜽) and their velocities, related to the mode shapes through the corresponding eigenvalue, so  𝜽𝑟= [𝜽𝑡 𝑟𝑠𝑟𝜽𝑡 𝑟]𝑡. This relationship can be expressed as a matrix as shown in Eq. (4), where 𝐒is the eigenvalue matrix, a diagonal matrix containing all the eigenvalues 𝑠𝑟and their complex conjugates 𝑠∗ 𝑟, Θis the eigenvector matrix whose columns are the eigenvectors  𝜽𝑟and their complex conjugates  𝜽∗ 𝑟ordered according to 𝐒, and Θis the modal shape matrix, whose columns are the mode shapes 𝜽𝑟and their complex conjugates 𝜽∗ 𝑟, also ordered according to 𝐒.  Θ=[Θ Θ𝐒](4) The diagonalization of the system in Eq. (2) can be performed thanks to the orthogonality properties with respect to the matrices 𝐀and 𝐁, expressed in Eq. (5). By taking advantage of the internal forms of 𝐀,𝐁and  Θ, the orthogonality properties can also be expressed in terms of the mass (𝐌), damping (𝐂), stiffness (𝐊) and mode shape (Θ) matrices, as shown in Eq. (6). Matrices  𝐀and  𝐁are both diagonal modal matrices, satisfying  𝐀−1  𝐁= −𝐒[6]. It is usual to scale the mode shapes so that matrix  𝐀is an identity matrix of dimension 2𝑛, and thus  𝐁equals −𝐒, a common assumption that is also adopted in this article.  Θ𝑡𝐀 Θ= 𝐀  Θ𝑡𝐁 Θ= 𝐁(5) Θ𝑡𝐂Θ+𝐒Θ𝑡𝐌Θ+Θ𝑡𝐌Θ𝐒= 𝐀 Θ𝑡𝐊Θ−𝐒Θ𝑡𝐌Θ𝐒= 𝐁(6) 3. Methodology The methodology proposed in this work, PhysEx, which is developed in the next subsections, provides a means to estimate physical models, and is composed of a set of three full-rank 𝑛×𝑛physical matrices, useful for the typical applications in structural engineering, from a set of 𝑚experimentally estimated modes. Note that the number of degrees of freedom of the resulting physical model 𝑛equals the number of monitored points during the measurement of the 𝑚modes, so the measured mode shapes are expressed only in those degrees of freedom. Both quantities usually differ, and this work focuses on the case in which 𝑚≤𝑛. The method is applied in two stages. In the first one, a system of linear equations is built to take into account the three main conditions that the sought matrices should meet. The first stage ends by imposing a nonlinear inequality condition, so an optimization process, which is the second stage, must be carried out to solve the problem. 3.1. Development of the system of equations The desirable set of physical matrices should fulfil certain conditions to provide an acceptable dynamic behaviour. In this work, the following three main conditions are considered: 1. The three physical matrices must be symmetric. This condition affects the internal form of the three matrices, forcing them to satisfy Eq. (7), where 𝐏may be any estimated physical matrix 𝐌𝑒,𝐂𝑒or 𝐊𝑒and subindices 𝑖and 𝑘stand for the row and column, respectively. 𝐏𝑖𝑘 −𝐏𝑘𝑖 = 0 ∀𝑖≠𝑘, (𝑖, 𝑘)≤𝑛(7) 2. The physical model must satisfy the eigenvalue problem in Eq. (3), which can be rewritten in terms of the physical matrices and the mode shapes as shown in Eq. (8). As the complex eigenvalues of a set of real-valued symmetric matrices always appear in conjugate pairs, this condition is enough to be met for the 𝑚pairs (𝑠𝑟,𝜽𝑟). (𝐌𝑒𝑠2 𝑟+𝐂𝑒𝑠𝑟+𝐊𝑒)𝜽𝑟= 0 𝑟≤𝑚(8) 3. If the previous conditions are met, the modal matrices  𝐀and  𝐁happen to be diagonal and satisfy  𝐀−1  𝐁= −𝐒. Therefore, the only remaining necessary condition is the mode shape scaling, i.e., the estimated physical matrices, together with the measured modal parameters, must satisfy Eq. (9). As in the previous condition, this one must only be fulfilled for the 𝑚pairs (𝑠𝑟,𝜽𝑟). 𝜽𝑡 𝑟𝐂𝑒𝜽𝑟+ 2𝑠𝑟𝜽𝑡 𝑟𝐌𝑒𝜽𝑟= 1 𝑟≤𝑚(9) Journal of Sound and Vibration 510 (2021) 116277 4 A. Magdaleno and A. Lorenzana Note that the positiveness of the matrices is not taken into account. This is further commented at the end of the first example. All the presented conditions, which are linear with respect to the physical properties, can be rewritten in the form of a system of equations, 𝐀𝑝𝐱𝑝=𝐛𝑝, where 𝐱𝑝contains the 3𝑛2elements of the physical matrices (the unknowns of the problem) and the internal form of 𝐀𝑝and 𝐛𝑝depends on the conditions being accounted for. After processing the symmetry conditions, the corresponding rows of 𝐀𝑝will contain 1and −1 in the appropriate locations; the associated terms in 𝐛𝑝will be 0. For the eigenproblem conditions, the corresponding terms in 𝐀𝑝will contain elements of different mode shapes, some of them multiplied by one eigenvalue; the corresponding terms in 𝐛𝑝will be 0 again. Finally, the terms of 𝐀𝑝associated to the scaling equations will be products of mode shape elements, some of them multiplied by their corresponding eigenvalue; the associated values in 𝐛𝑝are 1this time. Finally, note that after computing conditions 2 and 3, Eq. (8) and (9), many terms in the matrix 𝐀𝑝and vector 𝐛𝑝will be complex values, so, in general, the solution of the system will also be complex. This leads to complex physical matrices, which is something to be avoided to guarantee proper results when solving the system of differential equations in Eq. (1). To overcome this issue, the subset of complex equations 𝐀𝑝𝑐 𝐱𝑝=𝐛𝑝𝑐 can be split into their real and imaginary parts, so (Re (𝐀𝑝𝑐 )+ jIm (𝐀𝑝𝑐 ))𝐱𝑝= Re (𝐛𝑝𝑐 )+ jIm (𝐛𝑝𝑐 ), where the subindex 𝑐stands for the subset of complex-valued equations and Re () and Im () stand for real and imaginary parts, respectively. Given that 𝐱𝑝must be a real-valued vector, the previous complex equations must be satisfied for their real and imaginary parts separately, leading to Eq. (10), where all terms in 𝐀𝑞and 𝐛𝑞are real-valued and 𝐀𝑝𝑟 and 𝐛𝑝𝑟 account for the real-valued equations in 𝐀𝑝𝐱𝑝=𝐛𝑝. Note that this operation is to be carried out only with the eigenvalue and scaling conditions (that are inherently complex-valued), doubling their number of associated equations. Symmetry equations are not affected by this modification (and, thus, included in 𝐀𝑝𝑟 and 𝐛𝑝𝑟). 𝐀𝑞𝐱𝑝=𝐛𝑞 ⎡⎢⎢⎣ 𝐀𝑝𝑟 Re (𝐀𝑝𝑐 ) Im (𝐀𝑝𝑐 )⎤⎥⎥⎦ 𝐱𝑝=⎡⎢⎢⎣ 𝐛𝑝𝑟 Re (𝐛𝑝𝑐 ) Im (𝐛𝑝𝑐 )⎤⎥⎥⎦ (10) Summing up, the symmetry equations are associated to the out-of-diagonal terms of the three physical matrices, adding a total of 3𝑛(𝑛− 1)∕2 equations. There are 𝑛eigenproblem equations for each of the 𝑚modes, making a total of 2𝑛𝑚 eigenproblem equations, accounting for the modification in Eq. (10) as they are complex equations. Finally, there are 2 m real-valued scaling equations associated to the 𝑚diagonal elements of  𝐀. Therefore, the total number of necessary independent equations is, a priori, 3𝑛(𝑛−1)∕2+2𝑚(𝑛+1), which can only be exactly solved for the 3𝑛2unknowns only if the number of equations is less than or equal to the number of unknowns, leading to the restriction in Eq. (11), which couples the number of degrees of freedom 𝑛and the number of considered modes 𝑚. In any other case, the conditions can only be approximately met. 𝑚≤3 4𝑛(11) Note that the previous three conditions would also have been satisfied if the same approach had been directly applied to the whole set of orthogonality expressions in Eq. (6). However, in that case, the apparent number of equations would have been 8𝑚2 (twice (2𝑚)2), leading to a very different, and wrong, solvability condition to the relation between 𝑚and 𝑛. Finally, it is interesting to note that, for this methodology to be successfully applied, the required amount of modal information is limited by the solvability condition, and the usual belief that a better identification can be performed if the amount of information available is greater, cannot be applied here. If it is not met, unacceptable results may be obtained. So, in order to improve them, more DOFs need to be monitored to increase the order of the estimated physical matrices (𝑛) or the number of considered modes (𝑚) must be reduced, so that the estimated physical model will be valid in a narrower frequency band. 3.2. Solution to the system of equations A solution to the resulting system of equations can always be found if the restriction in Eq. (11) is met. However, if 𝑚is strictly lower than 3𝑛∕4, the system of equations is underdetermined and there is an infinite number of solutions, not all of them being suitable for engineering purposes. When dealing with incomplete models, one major problem that may arise is the characteristics of the exogenous modes, i.e., the modes that appear in the estimated model, which mathematically belong to it and can add unexpected and undesirable dynamic effects. This problem could reach the point of some exogenous modes having negative damping, which is terrible for time domain simulations. By controlling the effect of those exogenous modes, an optimal solution can be calculated. Let 𝐱𝑝𝑎 and 𝐱𝑝𝑏 be two particular solutions to the real-valued system of equations in Eq. (10), and 𝐱𝑝𝑑 their difference, so 𝐱𝑝𝑎 =𝐱𝑝𝑏 +𝐱𝑝𝑑 . Substituting 𝐱𝑝𝑎 in the system of equations leads to 𝐀𝑞𝐱𝑝𝑎 =𝐀𝑞(𝐱𝑝𝑏 +𝐱𝑝𝑑 ) = 𝐛𝑞and, given that 𝐱𝑝𝑏 also verifies it (𝐀𝑞𝐱𝑝𝑏 =𝐛𝑞), then 𝐀𝑞𝐱𝑝𝑑 = 0, which is the definition of the kernel of 𝐀𝑞. This means that, if the system of equations is underdetermined, so there is an infinite number of solutions, it is possible to obtain any solution to the system of equations from another one by adding a vector belonging to the kernel of 𝐀𝑞. A first solution can easily be obtained by means of the Moore–Penrose pseudo-inverse, which provides the solution with the minimum Euclidean norm. Then, the general solution can be expressed as shown in Eq. (12), where 𝐰𝑞contains a linear combination of the base vectors of ker (𝐀𝑞). 𝐱𝑝=𝐀+ 𝑞𝐛𝑞+ ker (𝐀𝑞)𝐰𝑞(12) The previous approach has a serious advantage in terms of the number of unknowns, having transformed a problem of 3𝑛2 variables into a problem of (3𝑛∕2 − 2𝑚)(𝑛+ 1) variables, the difference between the number of columns and the number of rows in Journal of Sound and Vibration 510 (2021) 116277 5 A. Magdaleno and A. Lorenzana 𝐀𝑞. Note that this value may be slightly higher if numerical linear dependencies appear in the system of equations. In any case, the number of terms in 𝐰𝑞is always fewer than 3𝑛2and always lead to a solution to the system of equations, i.e., a set of physical matrices containing, at least, the desired dynamic behaviour. This improvement can significantly reduce the computational effort and time when finding the final solution to the problem. Also note that, if the system of equations is overdetermined (because the condition in Eq. (11) is not met), then ker (𝐀𝑞)does not contain any vector and a solution of the system of equations can be found through 𝐱𝑝=𝐀+ 𝑞𝐛𝑞, which corresponds to the best fit to all the conditions in the least squares sense [23]. To get the final physical properties, an optimization problem is carried out. The objective is to obtain a solution of 𝐰𝑞that minimizes the interference of the exogenous modes, making sure that none of them have a negative damping value. Amongst the different objective functions that could be sought to achieve that, two of them are presented in the following section. 3.3. Objective functions 3.3.1. Exogenous natural frequency minimization It is known that if a physical model is directly obtained from the orthogonality conditions of an incomplete modal model [12,13], it contains exactly the desired modes with no exogenous modes, since their natural frequencies are equal to zero. This is due to the physical matrices being rank-deficient and, by definition, the last singular values of a rank-deficient matrix are zero, and so are the corresponding eigenvalues [23]. In that sense, the objective function proposed in this subsection is intended to minimize the natural frequencies of the exogenous modes in order to make the final physical model approach a rank-deficient one as much as possible. When a highly restrained optimization problem is solved, as it is the case here, it is unlikely that the objective function will be perfectly satisfied. For this reason, the odds of any exogenous natural frequency being exactly zero are really low and, in general, it can be assumed that the obtained physical matrices will not be rank-deficient. The exogenous natural frequencies may eventually remain close, or even inside, the frequency band of interest, affecting the dynamic behaviour of the model. To mitigate that, the natural frequencies are inversely weighted by their corresponding damping ratios to help dampen the exogenous eigenvalues and reduce their effect on the structural response. This objective function can be formulated as stated in Eq. (13), where the natural frequencies have been expressed as |𝑠𝑟|and the damping ratios as −Re (𝑠𝑟)∕|𝑠𝑟|, the |⋅|operator being the absolute value of a real number and the modulus of a complex number. 𝐽= 𝑛 ∑ 𝑟=𝑚+1 ||||||𝑠𝑟|2 −Re (𝑠𝑟)|||||(13) The outer absolute value is needed to prevent the negative and small Re (𝑠𝑟)values from becoming large negative 𝐽values, dominating the minimization process; although the negative parts of 𝑠𝑟are imposed as a nonlinear constraint in Section 3.4, due to the internal numerical constraint tolerances of the optimization algorithm. 3.3.2. Frequency response function approach Another approach is to minimize the error between two sets of Frequency Response Functions (or FRFs): the set calculated by means of the experimental modal model and the set associated to the estimated physical model. The FRF matrix 𝐇(𝜔)can be calculated from the experimental modal properties through  ΘΛ(𝜔) Θ𝑡[6], where Λ(𝜔)is a 2𝑚× 2𝑚 diagonal matrix containing the terms 1∕(j𝜔−𝑠𝑟), one for each eigenvalue including their complex conjugates, and 𝜔stands for the excitation frequency in rad/s. On the other hand, the relationship between the estimated physical matrices and the FRF matrix 𝐇𝑒(𝜔) is stated as (−𝐌𝑒𝜔2+ j𝐂𝑒𝜔+𝐊𝑒)−1, where a non-linear relationship can be appreciated. In the scenario in which the physical matrices exactly represent the same behaviour as the modal model in the frequency range of interest, both 𝐇(𝜔)and 𝐇𝑒(𝜔)should be equal for every 𝜔inside that frequency range. This means that 𝐇(𝜔)𝐇(𝜔)−1 𝑒=𝐈𝑛, where 𝐈𝑛is the identity matrix of dimension 𝑛. This fact leads to the expression in Eq. (14), which represents a linear relationship between the experimental modes and the physical matrices to be estimated. Note that this equation is not automatically met after Eq. (8) has been imposed due to the physical matrices 𝐌𝑒,𝐂𝑒and 𝐊𝑒also containing some exogenous modes, the influence of which (on the frequency range of interest) is to be reduced by means of the optimization procedure.  ΘΛ(𝜔) Θ𝑡(−𝐌𝑒𝜔2+ j𝐂𝑒𝜔+𝐊𝑒) = 𝐈𝑛(14) This relationship can be rewritten in the form of a system of linear equations 𝐀frf 𝐱𝑝=𝐛frf and, making use of Eq. (12), the cost function in Eq. (15) is obtained, where 𝑁stands for the number of elements in vector Γ. As it is a linear relationship, this expression is faster to compute than a regular nonlinear least squares optimization. It is also recommended to build 𝐀frf and 𝐛frf by selecting a sufficient number of well separated frequencies 𝜔inside the frequency range of interest, to make sure that an overdetermined system is built and that the exogenous modes do not affect any region of the frequency band of interest. 𝐽= 𝑁 ∑ 𝑖=1 |𝛤𝑖| Γ=𝐀frf ker (𝐀𝑝)𝐰𝑝+𝐀frf 𝐀+ 𝑝𝐛𝑝−𝐛frf (15) Journal of Sound and Vibration 510 (2021) 116277 6 A. Magdaleno and A. Lorenzana 3.4. The optimization strategy Although both previous objective functions can lead to exogenous natural frequencies not affecting the range of interest, it is still mathematically possible for their associated damping values to be negative. To avoid this, it is important to add a nonlinear constraint to the optimization process, forcing the algorithm to check the real parts of 𝑠𝑟on every iteration to make sure they remain negative. The addition of a nonlinear constraint may complicate the overall optimization process, since the algorithm must look for feasible points (i.e., vectors 𝐰𝑝that satisfy it) and minimize the objective function at the same time. Also, not providing a feasible initial point can worsen this issue, leading to bad optimization results when the number of variables is moderately large. To avoid this, the optimization procedure is split into two stages: 1. The optimization algorithm is first run including the nonlinear constraint, but with a null objective function (𝐽= 0). In this way, the objective function convergence is guaranteed from the beginning of the process and the algorithm can focus on finding a feasible point. This stage can be run several times using different randomly generated initial vectors to obtain a set of feasible initial points to be used in the second stage. 2. The second stage consists in running the optimization algorithm to minimize the selected objective function and starting from one of the feasible initial points found in the previous stage. With this approach, the algorithm can now focus on minimizing the objective function and just make sure it does not exit the feasible region from which it starts. The previous strategy does not ensure finding the global minimum, but provides different solutions that satisfy the constraints and minimizes the objective function inside the feasible region. Among them, the solution that provides the least objective function value is kept as the solution to the problem. Alternatively, a global optimization procedure could be applied, but this is not covered in this work. Here, the fmincon function is used to perform the different optimization. This function provides several gradientbased algorithms, such as the Interior-Point or the Trust-Region-Reflective, among others, from which the first is selected due to its versatility and ability to deal with medium to large sized optimization problems. All of them are prepared to deal with linear and nonlinear constraints, penalizing the objective function if any of them is not met. More details about the algorithm itself can be found in [24] and more details about the specific implementation of the methodology in MATLAB can be found in the provided code [22]. 3.5. Performance indicators In order to ease the comparisons between the different models that are to be obtained in Section 4, some indicators are proposed to account for the accomplishment of the various conditions commented in the previous subsections. 3.5.1. Symmetry indicator To measure the extent to which the estimated matrices are symmetric, the indicator in Eq. (16) is proposed, where ‖⋅‖𝐹stands for the Euclidean matrix norm (also known as the Frobenius norm) and the subindex 𝑑stands for the main diagonal of the matrix (e.g., 𝐌𝑒𝑑 = diag (𝐌𝑒)). This indicator accounts for the relative deviation of the upper triangular matrix with respect to the lower one. A smaller value of 𝜆𝑠𝑦 indicates a better symmetry and a value of 0 indicates a perfect symmetry. 𝜆𝑠𝑦 =1 3(1 2‖𝐌𝑒−𝐌𝑡 𝑒‖𝐹 ‖𝐌𝑒−𝐌𝑒𝑑 ‖𝐹 +1 2‖𝐂𝑒−𝐂𝑡 𝑒‖𝐹 ‖𝐂𝑒−𝐂𝑒𝑑 ‖𝐹 +1 2‖𝐊𝑒−𝐊𝑡 𝑒‖𝐹 ‖𝐊𝑒−𝐊𝑒𝑑 ‖𝐹)(16) It is important to note that, if one or more estimated physical matrices turn out to be diagonal, its corresponding denominator in Eq. (16) may become null, thus leading to numerical issues. As shown later, if the solvability condition in Eq. (11) is met, the terms of the estimated physical matrices will have no direct physical meaning and the strategy proposed in this work to obtain them will rarely lead to such a situation, unless it is enforced by some additional conditions. If such is case, the corresponding terms in Eq. (16) should be removed in order to avoid the aforementioned numerical issues. 3.5.2. Eigenvalues and eigenvectors indicator To measure how well the eigenvalues 𝑠𝑟, with 𝑟= 1..𝑚, are reproduced by the estimated physical model, the expression in Eq. (17) is proposed. This indicator calculates the mean relative distance between the expected complex eigenvalue 𝑠𝑟and the estimated one, 𝑠𝑟, which can be calculated as the closest one to 𝑠𝑟. When a perfect correspondence exists between the calculated and desired eigenvalues, then 𝜆𝑠= 0. 𝜆𝑠=1 𝑚 𝑚 ∑ 𝑟=1 |𝑠𝑟−𝑠𝑟| |𝑠𝑟|(17) The measurement of how well the mode shapes are replicated by the physical model can be done through the expression in Eq. (18), where the superindex 𝐻accounts for the conjugate transpose. This expression is based on the Modal Assurance Criterion (MAC) that is calculated for each pair of desired mode shape 𝜽𝑟and the corresponding calculated mode shape 𝜽𝑟. The first 1in Eq. (18) is included to make 𝜆𝜃vary between 0 (perfect correspondence) and 1 (no correspondence). 𝜆𝜃= 1 − 1 𝑚 𝑚 ∑ 𝑟=1 |𝜽𝐻 𝑟𝜽𝑟|2 (𝜽𝐻 𝑟𝜽𝑟)(𝜽𝐻 𝑟𝜽𝑟) (18) Journal of Sound and Vibration 510 (2021) 116277 7 A. Magdaleno and A. Lorenzana Fig. 1. Schema of the 3 DOF discrete model. Table 1 Natural frequencies (𝜔𝑟, in rad/s), damping ratios (𝜁𝑟, in %) and modal coordinates (𝜃𝑟𝑖, scaled to matrix 𝐴) of the original 3 DOF discrete model. 𝜔𝑟𝜁𝑟𝜃𝑟1𝜃𝑟2𝜃𝑟3 16.55 2.82 −0.025 − 0.810j −0.107 + 0.556j 0.111 − 0.103j 214.53 13.07 −0.427 − 0.171j −0.643 − 0.091j 0.558 + 0.235j 323.11 3.53 0.219 − 0.205j 0.047 − 0.489j 0.007 − 0.846j 3.5.3. Scaling indicator The last indicator accounts for the ability of the damping and mass matrices to provide a proper scaling of the modes, as stated in Eq. (6). The indicator in Eq. (19) is defined, which gets the value of zero only if the estimated matrices and the experimental mode shapes exactly satisfy Eq. (6). 𝜆𝑠𝑐 =1 𝑚 𝑚 ∑ 𝑟=1 |𝜽𝑡 𝑟𝐂𝑒𝜽𝑟+ 2𝑠𝑟𝜽𝑡 𝑟𝐌𝑒𝜽𝑟− 1|(19) 4. Application examples 4.1. 3 DOF discrete model The PhysEx procedure is now tested by means of the three-DOF discrete model shown in Fig. 1, with the following properties: 𝑚1= 2 kg, 𝑚2= 1.75 kg, 𝑚3= 1.5kg, 𝑐1= 0.3Ns/m, 𝑐2= 0.5Ns/m, 𝑐3= 4.1Ns/m, 𝑘1= 530 N/m, 𝑘2= 320 N/m and 𝑘3= 150 N/m. This ensemble of properties provide the physical matrices shown in Eq. (20) for a DOF vector 𝐪(𝑡) = [𝑥1(𝑡)𝑥2(𝑡)𝑥3(𝑡)]𝑡and the modal properties shown in Table 1. Note that the damping matrix is not proportional to the mass and the stiffness matrices, so the mode shapes are complex and general damping approaches must be used to estimate a physical model. For the frequency domain comparison, only accelerances, i.e., the relationship between a force input and an acceleration output, are computed assuming a force applied on the DOF 𝑥2. Time domain comparisons are made through the impulse response applying the input on the same DOF. 𝐌=⎡⎢⎢⎣ 2.00 0 0 0 1.75 0 0 0 1.50⎤⎥⎥⎦ 𝐂=⎡⎢⎢⎣ 0.80 −0.50 0 −0.50 4.60 −4.10 0 −4.10 4.10 ⎤⎥⎥⎦ 𝐊=⎡⎢⎢⎣ 850 −320 0 −320 470 −150 0 −150 150 ⎤⎥⎥⎦ (20) First, a complete model estimation is made. The estimated physical matrices, obtained using the information of Table 1, are shown in Eq. (21), where it is shown that they match well enough with the original ones. When the three modes are considered, the restriction in Eq. (11) concerning the number of modes is not met and the kernel of 𝐀𝑞has null dimension, so the system of equations is over determined and it can only be solved by pseudo-inverting matrix 𝐀𝑞. This implies that the provided solution is the best fit to all linear restrictions, only approximately satisfying them, which is evidenced in matrices 𝐌𝑒and 𝐂𝑒not being perfectly symmetric, as is also shown in the first row of Table 2. Here, the value of 𝜆𝑠𝑦 accounts for the lack of symmetry, the upper triangular matrices of the physical matrices being, on average, 6.38% different from the lower triangular ones. The rest of the indicators provide reasonably good results, but the indicator values are higher than the rest of the cases which do satisfy the Journal of Sound and Vibration 510 (2021) 116277 8 A. Magdaleno and A. Lorenzana Fig. 2. Accelerance comparison of the complete model for the DOFs (a) 𝑥1and (b) 𝑥3. restriction in Eq. (11). 𝐌𝑒=⎡⎢⎢⎣ 1.997 0.001 −0.002 0.002 1.745 −0.002 −0.003 −0.001 1.491 ⎤⎥⎥⎦ 𝐂𝑒=⎡⎢⎢⎣ 0.792 −0.478 −0.014 −0.477 4.692 −4.103 −0.012 −4.102 4.063 ⎤⎥⎥⎦ 𝐊𝑒=⎡⎢⎢⎣ 848.1 −318.9 −0.536 −318.9 −468.9 −149.3 −0.536 −149.3 149.0⎤⎥⎥⎦ (21) As stated by Lancaster [25], a complete and generally damped modal model only provides a symmetric and real-valued set of physical matrices if the mode shape matrix Θcan be expressed in terms of an orthogonal matrix 𝐆so Θ=Θ𝑅(𝐈𝑚− j𝐆), where the subscript 𝑅stands for the real part and 𝐈𝑚represents an identity matrix of dimension 𝑚×𝑚. In this case, the matrix 𝐆is shown in Eq. (22) and the product with its transpose is shown in Eq. (23), revealing that it is almost an orthogonal matrix, thus leading to almost symmetric and real-valued estimated physical matrices. However, structures not satisfying the solvability condition in Eq. (11) or the one proposed by Lancaster will not provide such meaningful results, as can be seen in the second example. 𝐆=⎡⎢⎢⎣ 0.992 0.122 0.008 −0.117 0.970 −0.215 −0.034 0.213 0.978 ⎤⎥⎥⎦ (22) 𝐆𝐆𝑡=⎡⎢⎢⎣ 0.999 0.001 0.000 0.001 1.001 0.001 0.000 0.001 1.003⎤⎥⎥⎦ (23) Fig. 2 shows the accelerance comparison between the original model and the estimated physical model in terms for the DOFs 𝑥1 and 𝑥3. As can be seen, there is a total correspondence between both responses in the frequency range of interest, so they can be said to be dynamically equivalent. This is a consequence of the good values obtained for the rest of the indicators in the first row of Table 2. Although the restriction in Eq. (11) is not met, in this case this is not so bad as the dynamics seems to be reproduced well enough. As will be shown in the next example, this is not usually true. The same conclusions can be drawn from the time domain comparison of Fig. 3, where the impulse response of DOF 𝑥3is compared. The intended application scenario of the proposed methodology is the one in which the number of considered modes satisfy Eq. (11), leading to incomplete models. In this case, PhysEx provides a set of three matrices that accurately satisfy the linear and nonlinear restrictions. However, the specific values that make up the matrices depend on the initial points used for the first stage of the optimization. This may lead to an infinite number of solutions, all suitable for engineering purposes, such as simulation or structural modification. Now, the estimated models resulting from considering only modes 1 and 2 are shown. Since two objective functions have been defined, two models are proposed: ‘‘1,2(f)’’ is associated to the FRF objective function and ‘‘1,2(w)’’ to the minimum exogenous frequency objective function. Eq. (24) shows the physical matrices associated to these models, where the subindex 𝑓stands for the Journal of Sound and Vibration 510 (2021) 116277 9 A. Magdaleno and A. Lorenzana Fig. 3. Time domain comparison of the complete model. Table 2 Values of the indicators calculated for each estimated physical model. Case 𝜆𝑠𝑦 𝜆𝑠𝜆𝜃𝜆𝑠𝑐 1,2,3 0.0638 3.04 × 10−6 6.34 × 10−7 2.16 × 10−3 1,2(f) 1.89 × 10−8 1.51 × 10−9 6.61 × 10−15 1.38 × 10−6 1,2(w) 4.53 × 10−10 1.51 × 10−9 6.55 × 10−15 1.38 × 10−6 1(f) 1.74 × 10−15 1.42 × 10−14 3.33 × 10−16 2.33 × 10−14 1(w) 2.26 × 10−16 3.45 × 10−15 0 1.00 × 10−15 model ‘‘1,2(f)’’ and the subindex 𝑤, for the model ‘‘1,2(w)’’. As can be seen, both sets of matrices are quite different from each other and from the original matrices of Eq. (20). It is interesting to note that the internal terms of these matrices do not have any direct physical meaning, much like a set of condensed matrices that lack a straightforward correspondence with the physical elements they represent. In addition, the estimated physical matrices have very different exogenous modes, which accentuates the differences in their internal terms. In spite of this, however, these physical matrices can still be treated as mass, stiffness and damping matrices, since they couple the accelerations, velocities and displacements of the considered DOFs, so they have a meaning as a whole. Finally, they are useful for the purposes of this work, since they contain the desired dynamic behaviour and, as will be shown later, they accurately predict the effect of mass modifications. The exogenous eigenvalues are −2.69 × 10−5 rad/s and −337.04 rad/s for ‘‘1,2(f)’’ and −6.05 × 10−5 rad/s and −1.41 × 10−4 rad/s for ‘‘1,2(w)’’. Their associated natural frequencies are well outside the frequency range of interest, although their damping ratio is 100% and their contribution to the time domain response would be low in any case. Note that two eigenvalues are obtained in each case because the size of the eigenvalue problem (Eq. (3)) is 2𝑛, so the two imposed eigenvalues determine a total of 4 (including their complex conjugates) with two exogenous ones remaining. Even if both models accurately represent the dynamic behaviour of interest, the FRF objective function lets the exogenous natural frequencies evolve towards large values or small values; whereas the minimum eigenvalue objective function imposes the evolution of the exogenous natural frequencies towards zero. In both cases, the four indicator values are very close to zero (rows 2 and 3 in Table 2), so they can be considered suitable models. 𝐌𝑒𝑓 =⎡⎢⎢⎣ 1.453 −0.787 0.830 −0.787 2.870 −0.592 0.830 −0.592 1.645 ⎤⎥⎥⎦ 𝐌𝑒𝑤 =⎡⎢⎢⎣ −92.90 57.19 −9.103 57.19 −32.74 5.506 −9.103 5.506 0.603 ⎤⎥⎥⎦ 𝐂𝑒𝑓 =⎡⎢⎢⎣ 383.4 −232.6 38.77 −232.6 142.3 −25.38 38.77 −25.38 6.434 ⎤⎥⎥⎦ 𝐂𝑒𝑤 =⎡⎢⎢⎣ −138.0 −41.95 61.71 −41.95 104.9 −53.14 61.71 −53.14 17.06 ⎤⎥⎥⎦ 𝐊𝑒𝑓 =⎡⎢⎢⎣ 592.2 104.1 −190.3 104.1 41.34 −0.412 −190.3 −0.412 108.5⎤⎥⎥⎦ 𝐊𝑒𝑤 =⎡⎢⎢⎣ 58.08 76.15 −32.31 76.15 102.1 −30.10 −32.31 −30.10 84.91 ⎤⎥⎥⎦ (24) Fig. 4 shows a comparison between the original accelerances and the ones corresponding to models ‘‘1,2(f)’’ and ‘‘1,2(w)’’ for the DOFs 𝑥1and 𝑥3. As expected, the estimated physical models accurately represent the dynamic behaviour of the retained modes. In the frequency domain, the curves associated to the estimated models approximate the original FRFs well enough in the surroundings of the first two modes, the ones under consideration. The same can be observed in the time domain comparison, Fig. 5, where the responses of the different estimated models are also very similar and the response of the original model is, to some extent, reproduced. Naturally, the transient contribution associated to the third mode is not present in the response of the estimated models, but once this contribution is damped in the response of the original model, all three behave in much the same way. Journal of Sound and Vibration 510 (2021) 116277 16 A. Magdaleno and A. Lorenzana Fig. 13. Accelerance comparison of the complete cantilever model, considering the modes 1 and 2, for the DOFs (a) 𝑣2and (b) 𝑣4. Fig. 14. Accelerance comparison of the modified incomplete models estimated by means of the eigenvalue approach considering the first two modes for the DOFs (a) 𝑥2and (b) 𝑥4. Table 7 Modified natural frequencies and damping ratios calculated for each estimated incomplete physical model compared to the modified original one. Case 𝜔1(rad/s) 𝜁1(%) 𝜔2(rad/s) 𝜁2(%) 𝜔3(rad/s) 𝜁3(%) Original 6.82 6.51 52.3 2.51 158.3 1.75 1,2,3(w) 6.82 6.62 52.4 2.78 161.0 2.13 1,2(w) 6.82 6.66 52.9 2.89 the conditions that the physical matrices should satisfy in order to provide a useful model and it has been demonstrated that, if the number of modes is less than 3∕4 the number of measured DOFs, the problem has infinite solutions; otherwise, only the best fit to all conditions can be obtained by pseudo-inverting the system matrix. To obtain a final solution for an incomplete model, an optimization process is carried out using different objective functions to show that both lead to solutions that satisfy the linear conditions, differing in the exogenous eigenvalues, although they always have an associated positive damping. The selection of the most suitable objective function must be done by the engineer for each particular case. Journal of Sound and Vibration 510 (2021) 116277 17 A. Magdaleno and A. Lorenzana The procedure was first applied to a discrete 3-DOF model, which has demonstrated the abilities of PhysEx to reproduce the original physical matrices when a complete model is considered, even if not all the linear constraints are perfectly satisfied. However, as expected, when an incomplete one is considered, both objective functions accurately reproduce the desired dynamic behaviour, but treating the exogenous eigenvalues differently. In those cases, the physical matrices differ, not only from the original ones, but also from each other. This is due to the presence of the exogenous modes, which are also different from one model to another, leading to the internal elements of the estimated physical matrices not having a direct physical meaning. However, these matrices can still be treated as mass, damping and stiffness matrices because they represent the mathematical coupling between the accelerations, velocities and displacements, respectively. In any case, they contain the desired modes and the frequency and time domain responses fit the original responses as far as they are able to. As a consequence, the same mass modification can be applied to both the original and the estimated model and very similar effects are predicted from each one. Finally, a discretized cantilever beam has been analysed. Apart from the conclusions that have also been drawn from the previous example, this example has shown the differences between the dynamic condensation procedure, which does not exactly respect the original modes but provides a set of perfectly symmetric matrices, and PhysEx, which always tries to simultaneously meet all the linear constraints, even when the solvability condition is not met, resulting in different physical models. This example has also shown a situation in which the number of modes is exactly 3∕4 the number of DOFs. It may happen that, even in those cases, some numerical linear dependencies lead to an underdetermined system of equations and the optimization procedure needs to be applied, although for a small number of variables. This fact could potentially lead to situations in which the solvability condition is not theoretically met but, thanks to the numerical characteristics of the problem, the system ends up being underdetermined. CRediT authorship contribution statement Alvaro Magdaleno: Formal analysis, Software, Validation, Investigation, Resources, Writing - original draft. Antolin Lorenzana: Conceptualization, Methodology, Writing - review & editing, Supervision, Funding acquisition. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors wish to acknowledge the Spainsh State Research Agency, Ministerio de Economía y Competitividad, Spanish Government, for the partial support through the RTI2018-098425 Research Project and the Ministerio de Educación, Cultura y Deporte, Spanish Government, for the FPU16/01339 predoctoral grant. Reproducibility and open source policy The authors of this work have a consistent policy of making science and engineering codes available openly, in the interest of reproducibility. The entire code and data that was used to obtain the present results is available online from [22] and usage is licenced under CC-BY.4.0. Appendix. Model matrices This appendix contains some physical matrices, not necessary to follow the main ideas of this article. Journal of Sound and Vibration 510 (2021) 116277 18 A. Magdaleno and A. Lorenzana A.1. 8 DOFs original model of the cantilever beam 𝐌= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0.216 0.037 −0.009 0.006 0.009 −0.002 0.037 0.009 0.216 0.037 −0.009 −0.009 −0.002 0.006 0.009 −0.002 0.037 0.009 0.216 0.037 −0.009 −0.009 −0.002 0.006 0.009 −0.002 0.037 0.009 0.108 −0.015 −0.009 −0.002 −0.015 0.003 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ 𝐂= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0.320 0.320 0.320 0.320 0.320 0.320 0.320 0.320 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ 𝐊= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 16128 −8064 4032 5376 −4032 1344 −8064 −4032 16128 −8064 4032 4032 1344 5376 −4032 1344 − 8064 −4032 16128 −8064 4032 4032 1344 5376 −4032 1344 −8064 −4032 8064 −4032 4032 1344 −4032 2688 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (A.1) A.2. Matrices of the estimated cantilever model considering 3 modes Eq. (A.2) shows the matrices obtained by applying PhysEx with the FRF objective function and Eq. (A.3) shows the resulting matrices after applying PhysEx with the minimum eigenvalue objective function. 𝐌𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 0.091 0.131 −0.016 0.025 0.131 0.225 −0.031 −0.017 −0.016 −0.031 0.367 0.026 0.025 −0.017 0.026 0.085 ⎤⎥⎥⎥⎥⎦ 𝐂𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 1.451 −1.564 1.135 −0.356 −1.564 1.686 −1.224 0.384 1.135 −1.224 0.889 −0.279 −0.356 0.384 −0.279 0.088 ⎤⎥⎥⎥⎥⎦ × 103 𝐊𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 1.270 −0.530 −0.110 0.130 −0.530 0.400 −0.010 −0.078 −0.110 −0.010 0.047 −0.014 0.130 −0.078 −0.014 0.024 ⎤⎥⎥⎥⎥⎦ × 103 (A.2) 𝐌𝑒𝑤 =⎡⎢⎢⎢⎢⎣ −367.7 396.5 −287.7 90.34 396.5 −427.0 8.378 −97.36 −287.7 310.1 −2.696 70.67 90.34 −97.36 70.67 −22.09 ⎤⎥⎥⎥⎥⎦ 𝐂𝑒𝑤 =⎡⎢⎢⎢⎢⎣ −2.758 1.671 −0.307 −0.098 1.671 −0.396 −0.688 0.426 −0.307 −0.688 1.209 −0.532 −0.098 0.426 −0.532 0.215 ⎤⎥⎥⎥⎥⎦ × 103 𝐊𝑒𝑤 =⎡⎢⎢⎢⎢⎣ 0.321 0.586 −0.185 −0.103 0.586 1.081 −0.359 −0.191 −0.185 −0.359 0.149 0.072 −0.103 −0.191 0.072 0.037 ⎤⎥⎥⎥⎥⎦ × 103 (A.3) Journal of Sound and Vibration 510 (2021) 116277 19 A. Magdaleno and A. Lorenzana A.3. Matrices of the estimated cantilever model considering 2 modes Eq. (A.4) shows the matrices obtained by applying PhysEx with the FRF objective function and Eq. (A.5) shows the resulting matrices after applying PhysEx with the minimum eigenvalue objective function. 𝐌𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 5.564 −6.016 4.680 −1.638 −6.016 5.799 −3.536 1.090 4.680 −3.536 2.038 −0.386 −1.638 1.090 −0.386 0.152 ⎤⎥⎥⎥⎥⎦ 𝐂𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 333.7 −498.5 444.6 −156.0 −498.5 730.2 −645.1 224.9 444.6 −645.1 569.1 −198.1 −156.0 224.9 −198.1 69.46 ⎤⎥⎥⎥⎥⎦ 𝐊𝑒𝑓 =⎡⎢⎢⎢⎢⎣ 758.9 578.7 100.5 −340.9 578.7 420.7 −166.5 −76.57 100.5 −166.5 325.3 −140.7 −340.9 −76.57 −140.7 162.8 ⎤⎥⎥⎥⎥⎦ (A.4) 𝐌𝑒𝑤 =⎡⎢⎢⎢⎢⎣ −55.68 −55.15 106.9 −46.54 −55.15 15.71 15.29 −9.912 106.9 15.29 −89.36 43.37 −46.54 −9.912 43.37 −20.51 ⎤⎥⎥⎥⎥⎦ 𝐂𝑒𝑤 =⎡⎢⎢⎢⎢⎣ −111.2 −97.41 75.98 −5.363 −97.41 −62.22 108.2 −40.36 75.98 108.2 −59.77 −5.257 −5.363 −40.36 −5.257 18.44 ⎤⎥⎥⎥⎥⎦ × 103 𝐊𝑒𝑤 =⎡⎢⎢⎢⎢⎣ 298.0 395.9 21.01 −177.1 395.9 530.0 38.28 −231.8 21.01 38.28 18.47 −0.757 −177.1 −231.8 −0.757 107.5 ⎤⎥⎥⎥⎥⎦ × 103 (A.5) References [1] S. Lammens, W. Heylen, P. Sas, Model updating and structural optimisation of a tennis racket, Mech. Syst. Signal Process. 6 (5) (1992) 461–468, http://dx.doi.org/10.1016/0888-3270(92)90069-U. [2] E. Simoen, G. De Roeck, G. Lombaert, Dealing with uncertainty in model updating for damage assessment: A review, Mech. Syst. Signal Process. 56–57 (2015) 123–149, http://dx.doi.org/10.1016/j.ymssp.2014.11.001. [3] N.A.Z. Abdullah, M.S.M. Sani, M.M. Rahman, I. Zaman, A review on model updating in structural dynamics, IOP Conf. Ser. Mater. Sci. Eng. 100 (2015) 012015, http://dx.doi.org/10.1088/1757-899X/100/1/012015. [4] M.I. Friswell, J.E. Mottershead, Finite Element Model Updating in Structural Dynamics, in: Solid Mechanics and its Applications, vol. 38, Springer Netherlands, Dordrecht, 1995, http://dx.doi.org/10.1007/978-94-015-8508-8. [5] V. Arora, Comparative study of finite element model updating methods, J. Vib. Control 17 (13) (2011) 2023–2039, http://dx.doi.org/10.1177/ 1077546310395967. [6] D.J. Ewins, Modal Testing: Theory, Practice, and Application, second ed., Research Studies Press, 2000. [7] N. Maia, J. Silva, Theoretical and Experimental Modal Analysis, in: Engineering Dynamics Series, Research Studies Press, 1997. [8] M. Ouisse, E. Foltête, Identification of Reduced Models from Optimal complex eigenvectors in Structural Dynamics and Vibroacoustics, in: Vib. Struct. Acoust. Anal., Springer Netherlands, Dordrecht, 2011, pp. 303–327, http://dx.doi.org/10.1007/978-94-007-1703-9_11. [9] S. Mukhopadhyay, H. Lus, R. Betti, Structural identification with incomplete instrumentation and global identifiability requirements under base excitation, Struct. Control Health Monit. 22 (7) (2015) 1024–1047, http://dx.doi.org/10.1002/stc.1732. [10] A. Berman, W.G. Flannelly, Theory of incomplete models of dynamic structures, AIAA J. 9 (8) (1971) 1481–1487, http://dx.doi.org/10.2514/3.49950. [11] A. Berman, System identification of structural dynamic models Theoretical and practical bounds, in: 25th Struct. Struct. Dyn. Mater. Conf., American Institute of Aeronautics and Astronautics, Reston, Virigina, 1984, http://dx.doi.org/10.2514/6.1984-929. [12] M. Link, Theory of a method for identifying incomplete system matrices from vibration test data, Z. Flugwiss. Weltraumforsch. 9 (2) (1985) 76–82. [13] H.-P. Chen, N. Bicanic, Assessment of damage in continuum structures based on incomplete modal information, Comput. Struct. 74 (5) (2000) 559–570, http://dx.doi.org/10.1016/S0045-7949(99)00062-0. [14] S.R. Ibrahim, Dynamic modeling of structures from measured complex modes, AIAA J. 21 (6) (1982) 898–901, http://dx.doi.org/10.2514/3.8168. [15] C.-P. Fritzen, Identification of mass, damping, and stiffness matrices of mechanical systems, J. Vib. Acoust. Stress Reliab. Des. 108 (1986) 9–16, http://dx.doi.org/10.1115/1.3269310. [16] K.-S. Kim, Y.J. Kang, J. Yoo, Structural parameters identification using improved normal frequency response function method, Mech. Syst. Signal Process. 22 (8) (2008) 1858–1868, http://dx.doi.org/10.1016/j.ymssp.2008.02.001. [17] T. Kasai, M. Link, Identification of non-proportional modal damping matrix and real normal modes, Mech. Syst. Signal Process. 16 (6) (2002) 921–934, http://dx.doi.org/10.1006/mssp.2001.1478. [18] A. Srikantha Phani, J. Woodhouse, Viscous damping identification in linear vibration, J. Sound Vib. 303 (3–5) (2007) 475–500, http://dx.doi.org/10. 1016/j.jsv.2006.12.031. Journal of Sound and Vibration 510 (2021) 116277 20 A. Magdaleno and A. Lorenzana [19] S. Adhikari, A.S. Phani, Experimental identification of generalized proportional viscous damping matrix, J. Vib. Acoust. 131 (1) (2009) 011008, http://dx.doi.org/10.1115/1.2980400. [20] M. Prandina, J.E. Mottershead, E. Bonisoli, An assessment of damping identification methods, J. Sound Vib. 323 (3–5) (2009) 662–676, http://dx.doi.org/ 10.1016/j.jsv.2009.01.022. [21] A. Bajrić, J. Høgsberg, Identification of damping and complex modes in structural vibrations, J. Sound Vib. 431 (2018) 367–389, http://dx.doi.org/10. 1016/j.jsv.2018.05.048. [22] A. Magdaleno, A. Lorenzana, Physex: A novel procedure to estimate full-rank physical matrices of a structure from an incomplete modal model, Mendeley Data (2020) http://dx.doi.org/10.17632/r9pbp4zvwm.1. [23] S.L. Campbell, C.D. Meyer, Jr., Generalized Inverses of Linear Transformations, Dover Pub. Inc., New York, 1991, p. 272. [24] R.H. Byrd, M.E. Hribar, J. Nocedal, An interior point algorithm for large-scale nonlinear programming, SIAM J. Optim. 9 (4) (1999) 877–900, http://dx.doi.org/10.1137/S1052623497325107. [25] P. Lancaster, U. Prells, Inverse problems for damped vibrating systems, J. Sound Vib. (2005) http://dx.doi.org/10.1016/j.jsv.2004.05.003. [26] Z.-Q. Qu, Model Order Reduction Techniques with Applications in Finite Element Analysis, Springer London, London, 2004, p. 369, http://dx.doi.org/10. 1007/978-1-4471-3827-3.