A novel optimization framework using frequency-based substructuring for estimation of linear bolted joint stiffness and damping

Estimating the stiffness and damping coming from a joint is a major challenge. This work proposes a novel framework to estimate the unknown parameters using frequency-based substructuring and a maximum a posteriori optimization approach. The idea is to minimize the discrepancy between the measured responses of an assembled system and the responses obtained by coupling the measured responses of the individual substructures with a joint model. The system’s dynamics are assumed to be linear. However, the FRFs of the coupled system (which drive the objective function) change nonlinearly with the linear joint parameters, making it a nonlinear optimization problem. Therefore, the joint parameters are estimated iteratively using the Gauss minimization method, a well-established linearization scheme. The framework is numerically validated on simple mdof systems


Introduction
Determining and predicting joint stiffness and damping is a challenge [1].Over the last decades, extensive studies and innovative approaches using frequency-based methods have been proposed [2][3][4][5][6][7][8][9][10].Despite the significant progress, improvements are still needed to robustly identify and estimate joint parameters when dealing with complex industrial machinery and structures.In many applications, it would undoubtedly be a great help in a design phase if predictions of stiffness and damping of a bolted connection were available or if it was possible, with a few simple experiments, to predict the stiffness and damping of a given structure.As the literature shows [11,12], considering complex structures as a set of simpler substructures is a good starting point when identifying a joint.Accurate models exist for substructures with well-defined boundary conditions.Frequency-Based Substructuring (FBS) [13,14] makes it straightforward to rigidly connect well-defined structures in the frequency domain and predict the assembled response.However, when real parts get assembled, it is more complex.A real connection is rarely rigid.It will be flexible, but the question is, how flexible?The mass of a bolt and nut will also affect the response and cannot readily be neglected.The primary challenge with Table 1 *More methods and variations exists.This selection provides a background for the different choices to be made when estimating joint parameters.interface refers to measurements near the joint interface, and internal refers to measurements away from it.a In [23] an unknown cable is estimated and not a joint, but dynamic decoupling is applied.[16] ✓ ✓ ✓ Updating  J () FBS dynamic decoupling (algebraic) [10] ✓ ✓ ✓ Direct  J () FBS dynamic decoupling (negative substructuring) [23] a Inverse substructuring [18] ✓ Direct  J () all joint estimations is that the joint or connection only exists when two structures are assembled, so it cannot be isolated physically and measured alone.This work aims to find the best joint stiffness and damping parameters that, when coupling two substructures in the frequency domain, result in frequency responses as close as possible to those obtained when measuring the real assembled structure.The principal novelty in this work is the combination of a modified nonlinear maximum a posteriori (MAP) estimator with FBS, while applying the virtual point transformation scheme (VPT) [11,15] to make the method applicable for any given structure with a single bolted connection.
As aforementioned, extensive research on this topic exists already.Generally, it is important to distinguish between joint identification and joint parametrization.State-of-the-art methods, such as those proposed in [10,16], identify the joint admittances, leaving it as a subsequent optional step to estimate the joint parameters.In [6], the joint damping and stiffness are estimated per frequency line and subsequently averaged using statistical methods.In this work, it is done differently.An underlying joint model is chosen beforehand as an initial condition, and the optimization is based on that model.With an underlying joint model, it is possible to find the parameters giving the best correlation between the assembled system's measurements and the independent substructures coupled with that joint model.A benefit of parametrizing before optimizing is that it works as a regularization: it restricts the optimization process to this particular model, which, in case of poor data, helps to avoid divergence.The approach shares traits with several state-of-the-art methods, e.g., the method by S. Tol and H.N. Özgüven [10] (with similarities to [2] by J.S. Tsai et al.), and with the SEMM-based (System Equivalent Model Mixing) optimization method developed by S.W. Klaassen et al. [16,17].All require the same set of measurements.In the method proposed here, and in [10,16], the joint admittances are found without measuring the interface DoFs of the assembled structure, as opposed to the quasi-static decoupling method, known as inverse substructuring [18][19][20][21][22], that requires measurements near the interface on the assembled structure.By inverse substructuring we refer to methods needing only measurements of the assembled substructure, but requiring a quasi-static interface.
Table 1 provides a brief overview of the needed measurements and the obtained output for a selection of state-of-the-art methods in the hope that it clarifies where this novel MAP-FBS method is similar and distinctive compared to existing procedures.It is worth noting that all the joint estimation techniques represent compromises.The dream scenario would naturally be to require a minimum of measurement data and avoid complex iterative procedures, but at the same time, in the end, obtain parametrized values of the joint.Having it all is impossible, but this novel method is a compromise in which the interface DoFs of the assembled structure are unnecessary (as in SEMM and algebraic FBS decoupling).However, unlike those, this method directly parametrizes the joint, so no subsequent steps are necessary.The ''Direct'' methods in the table use no optimization to get an admittance (or impedance) as a function of the frequency of the joint.However, extending those methods with optimization is probably possible.By negative substructuring, we refer to the classical FBS coupling equation, using the negative admittance of a substructure to decouple it.Theoretically, FBS decoupling (algebraic and negative substructuring) works well, but in the presence of noise, it encounters problems due to the inversion of matrices.To try to overcome this problem, Wang [24] considered a method for choosing the best data, proposing an algorithm to reduce error when the noise is non-Gaussian or if the dynamic stiffness of the joint is much bigger than the one of the substructures it is attached to.The method is based on the FBS formulation from [25].Naturally, the effectiveness of all identification methods will depend on the quality of measurements and what is done to try to improve them.
Most methods and studies focus on joint stiffness estimation, as estimating damping is more challenging.In [10], the authors determine damping per frequency line based on real experimental data.The estimates vary by several hundred percentages over this range, but the compromise value they find provides a good agreement with reference data.In [8], joint damping is attempted to be estimated from experimental data.However, the estimates are mostly non-physical (e.g., negative damping values) due to an ill-conditioned system and a least-square approach that does not consider the physics.
The accuracy of each method in Table 1 will depend on the nature of the problem.It cannot be generally stated which method is better.If, for example, it is difficult to place sensors near the interface on the assembled structure, inverse substructuring and classical FBS decoupling are not good choices, but if the substructures, on the contrary, cannot easily be measured alone, inverse substructuring would be a good choice.The different methods are also not one-to-one comparable, as either different measurement data is needed, or the output of the methods differs.Methods yielding joint admittance  J require an additional step to derive joint model parameters.Here, signal filtering plays a critical role, and the quality of results depends on the user's proficiency in signal processing.A proper quantitative comparison requires substantial work to treat all methods fairly.Various structures and applications should be tested to reveal the most efficient method for each case.While a broad comparative study can help determine if one method performs better overall, this study's focus lies in an earlier stage, deriving and validating the novel MAP-FBS method.Therefore, a comprehensive comparison is better suited for a separate study.However, regardless of the chosen method, several general features drive the final quality in FBS-based joint estimation.These can be categorized into: • Measurement noise • Measurement uncertainty (positioning of sensors and excitation) • Visibility in measurements: Is the joint active in the measurements?
• The representativeness of the joint model (the consequence of simplified assumptions) • Strength of nonlinearity: Is it reasonable to neglect nonlinearity?• Joint stiffness/damping magnitude in comparison to other sources of stiffness/damping Measurement noise cannot be avoided, but options exist to improve and lessen their effect.Measurement uncertainty also follows with any experimental campaign.Hammer impacting in an exact predefined position and direction is impossible.A sensor's attachment position cannot be controlled entirely, and when using accelerometers for sensors, the exact measuring point within the attachment surface is also unknown.These uncertainties must be accepted but can be better understood when conducting uncertainty studies.E.g., by slightly varying the defined positions and directions and observing the effect on the joint parameters.
Visibility refers to the effect of the chosen measurement positions and the measured frequency range: Even the best method cannot estimate inactive parameters in the data.When planning an experiment, it can be helpful, via a FEM tool, to identify, with the use of a spring model for the joint, which modes are sensitive to a change in joint parameters.The measurement campaign should be designed so that the effect of the joint is visible in the recorded data.The fourth bullet point refers to the applied joint model.Any joint model is naturally a simplification.On the other hand, assuming a too complex model might lead to unreliable estimates and over-fitting.Using the VPT assumption, as in this work, has the benefit that no matter the structure, one can approximate the joint parameters as a 6-DoF connection.The VP makes it possible to compare the joint stiffness of entirely different structures.On the other hand, the VP comes with limitations.The flexible 6-DoF VP connection entirely ignores the specific geometry of the bolted connection, and does not model the fine details of the real contact surface.It is a mathematical abstraction, for which we can represent any bolted connection via translational and rotational linear stiffness (and damping) in all spatial directions.
In this work, the joint stiffness and damping are assumed linear, which is a fair assumption for stiff joints.However, one should always consider if the nonlinearity is so strong that a linear model will not represent the system sufficiently.Measured FRFs of an assembled structure can quickly reveal nonlinearity.If these have good coherences and clear straight resonance peaks, the nonlinearity from the joint is negligible.In case of large vibration amplitude, it is recommended to proceed with caution, as that in many instances induces non-negligible nonlinear behaviour in the joint due to changes in contact forces [26].
The last issue is related to magnitudes.Suppose the joint damping is minimal compared to other sources (e.g., damping from the elastic bands from which the structure is hanging).In that case, estimating accurate values might be difficult (and irrelevant).
The manuscript is organized as follows: first, a brief overview of frequency-based substructuring covering two different joint models and coupling approaches, along with an introduction of the Virtual Point, a mathematical abstraction allowing coupling of complex 3D structures.Next follows the derivation of the novel MAP-FBS method in Section 3. In Section 4, the novel method is validated using two numerical examples, and the effect of noise on the algorithm is discussed.The method is tested on real experimental data, using the benchmark structure known as the AM structure, to check the application potential.The results of the experimental investigation are presented and discussed in Section 5. Section 6 provides the concluding remarks.

Frequency-based substructuring and joint modelling
Numerous publications describe the method of frequency-based substructuring (FBS) for coupling dynamic systems [14,15].This section provides a short recap to motivate the optimization procedure presented in the next section.The first step is to combine the two systems' frequency response functions in a block matrix form.By defining where the two systems are connected, the response of the coupled system can be obtained.In terms of equations, the starting point is: where  is the impedance matrix,  the displacement vector,  a signed Boolean matrix mapping the matching DoFs and  = null(),  the external forces and  the interface forces.
This work applies the method known as dual coupling using the Lagrange Multiplier (LM) and considers three different connection interpretations, between two given substructures: Case 0, they are rigidly connected (this is just for reference purposes), Case 1, they are connected by springs and dampers, and Case 2, the connection is a substructure itself, with springs, dampers and masses [27].Starting with the most straightforward formulation, Case 0, and applying the dual assembly approach leads to:  where the block matrix  A|B holds the admittances of both substructure A and B. Note, the admittance functions are always function of , but the () is left out for compactness.By inserting the second line of Eq. ( 2) into the first and reorganizing, the FRFs' of the coupled response can be obtained as [15]: Note, ỸAB contains twice the connection DoFs.Subsequently the redundant DoFs may be removed if wished.Now, introducing flexibility (and potentially damping) in the connection, thus moving on to Case 1, the interface compatibility (second equation in (2)) is now replaced by a interface admittance  =  j .Eliminating again the interface forces  from the equation leads to: where  j is the flexibility matrix at the DoFs of the connection.This formulation is a quasi-static coupling and is often referred to as the weakened compatibility formulation.Fig. 1 illustrates the simplest interpretation of Case 1.For easy visualization, the example is one-dimensional, but (4) is applicable to multi-directional structures.For the next steps, with complex structures, it is useful to denote the DoFs connecting the two structures as interface DoFs and the remaining as internal DoFs.We return to how to handle 3D structures in the following Section 2.1.Note, with this quasi-static LM approach, it is not possible to have cross-coupling terms, as the interpretation of the interface forces  =  T  is only satisfied if each DoF only has one other contact DoF.
In the last case, Case 2, the joint is interpreted as a separate substructure, which requires an extension of the block matrix in (2): where  J is the admittance of the separate joint substructure.The capital J indicates the joint is a separate substructure.Fig. 2 shows the simplest interpretation of Case 2. The coupling follows the same form as in (3), but  is new, matching the coupling DoFs of this setup, where the block matrix includes now all three substructures: In the lab, it is, in most cases, easy to measure internal DoFs of an assembled setup  AJB mea , as well as the responses of the individual substructures  A and  B .However, measuring  J is normally impossible.It does not exist on its own.This work aims to use relation (4)/(6) (depending on the case) with a nonlinear optimization approach to find the best model parameters for  J (or  j ).E.g., assuming the left-hand side of ( 6) is a known quantity, thus replacing ỸAJB with the actual measured response  AJB mea , and everything on the right-hand side, except  J , is also known, then with nonlinear regression it is possible to crunch the numbers backwards to find the parameters for  J that fulfil the FBS equation the best.3. VP substructure interpretations.Despite the one-dimensional visualization, the VP is a three-dimensional abstraction, providing translational and rotational stiffness in all spatial directions.Note, for every spring there could also be a damper.The superscripts A and B indicates the DoFs for substructure A and B, respectively.Fig. 4. A real structure interpretation.For simplicity here illustrated in 2D, but for an actual experiment, the structure is in 3D and the impact points and sensors are placed along all three dimensions.As for the mDoF system, the impacts and sensors on Structure A and B away from the interface are named the internal DoFs.Those near the interface are the interface DoFs.The interface DoFs are transformed into the virtual point with 6-DoFs.

Virtual point transformation for analysing complex 3D structures
When handling complex structures, it is challenging to enforce collocated measurements at each DoF of two connecting subsystems.Another way of coupling complex interfaces is to use a virtual point (VP) [15].The idea is to connect the structure via a single virtual point only.For a 3D structure, a virtual point simplifies the connection of two structures to a 6-DoF connection, denoted : three translations (  ,   ,   ) and three rotations (   ,    ,    ).As it is impossible to measure the admittances directly at the VP, we measure near it, and by simple geometrical conditions transform the measured accelerations to get the VP accelerations.That can only be done by assuming that there is no flexibility between the measurement point and the VP.The standard procedure is to measure at several locations near the joint, obtaining an overdetermined system.The transformation is then the least square solution to all the measurement points.A VP assumes a rigid relation between the VP and the sensors near the interface, which means a VP interpretation is not valid for higher frequencies.The frequency range up to which the assumption of rigidity of the interface that is underlying the VP is valid can be estimated, for instance, by checking the consistency of the measurements [28] or by comparing the dynamics obtained with a FEM model for the true interface and for the interface made rigid.Cf. [15] for more details on the VPT.
The simplest option is to have a virtual point to rigidly connect the two structures, Case 0. It is also possible to model the connection so that the two substructures are connected via the virtual point by springs and dampers, as Case 1 in Fig. 3 illustrates.Lastly, the two substructures can also be connected via the virtual point with an independent substructure, with mass/inertia, as well as stiffness, as Case 2 in Fig. 3 demonstrates.As mentioned in the Introduction, the six DoFs VP is a 3D mathematical abstraction of the actual connection.The micro-scale contact surface details are not taken into account.The connection for a real structure is always defined as a VP, with 6 × 6 DoFs, but the number of internal DoFs depend on the system.Often is it easier to make many impacts, and only have a limited number of sensors.The number of measurements on A and B is also not necessarily the same.
The principle for coupling is the same as for the simpler mDoF systems shown in the previous section.Fig. 4 shows which measurements are needed.It is necessary to excite and measure a selected number of internal DoFs on both A and B, as well as to place a number of sensors and make some impacts near the interface.The measurement points close to the connection are denoted interface DoFs and are used to perform the virtual point transformation.After finding the translations and rotations in the virtual point, using simple geometrical relations [15], the connection between A and B will only be described via the VP.Then, the connection can be interpreted exactly as Case 0, 1, and 2 for the simple systems.

The maximum a posteriori likelihood (MAP) adapted to frequency-based substructuring
An objective function to be minimized is set up, inspired by the MAP-estimator [29]: where  are the unknown joint parameters, p is the a posteriori estimate (thus the initial guess of parameters),  AJB mea is a vectorization of the measured frequency-responses  AJB mea of the assembled structure at the internal DoFs.The order of how the responses  AJB mea are transformed into a vector is not important.Naturally, it has to be consistent, so the estimated responses ŷAJB () and measured responses  AJB mea are given in column format in the same way.The first term expresses a nonlinear least-square problem for the joint parameter , describing the difference between estimated and measured responses.The second term (called the a posteriori term since it uses an estimate p of the parameters) works as a regularization to avoid obtaining estimated responses that are non-physical or unrealistic joint stiffness and damping.E.g., negative stiffness. and  are weight matrices (symmetric and positive definite); ideally, for a proper MAP, this is the covariance of the measurements and the parameters, respectively.However, as those values are rarely available, the expected variance is used instead.
Note, a slight adjustment is introduced, compared to the classical MAP formulation [29]: As measured frequency response functions are complex functions, the first term minimizes the difference between the measured and estimated responses, but by taking the adjoint of the difference.The adjoint operation, denoted ()  is the complex conjugate of a transposed vector, ((⋅) T ) * .The first term will then be real and match the second, also real, term.
Minimizing the objective function in ( 7) is a nonlinear regression problem, as the dependent variable ŷAJB in ( 7) is nonlinear in the unknown parameters .This introduces significant complexity compared to linear regression.Although the frequency responses are assumed linear and the substructuring coupling is linear, the unknown parameters describing joint stiffness are nonlinearly related to the FRFs in Eqs. ( 4)/(6) (Case 1/Case 2).The objective function is therefore deliberately chosen; the minimizing solution would diverge or end up at local minima representing non-physical values without the second term in (7) due to this nonlinearity in terms of parameters.
The objective function is a choice among many available options.Slight variations were tested, such as using the logarithm of the responses, but first evaluations showed that it did not change the effectiveness significantly.In future work, variations may be interesting to study.This objective function can, as presented here, in principle, be implemented as per frequency line optimization or over all the measured frequencies.The only formal difference is the length of the vector  AJB mea .If the optimization is made per frequency line, the estimated joint parameters will be a function of frequency.Computation-wise, this would be costly, as one has to go through the iterative process for all frequencies instead of just one time.Otherwise, the cost function could possibly be adapted to apply a pre-defined function of frequency (()) to optimize over the entire frequency range in one step.In this work, we assume the joint parameters to be constant over the frequency range, and the optimization is global, based on all frequency lines at once.Fig. 5 illustrates what DoFs are used in the cost function  when we consider a real structure.The response ŷAJB mea corresponds only the internal DoFs, not those at the interface, as we assume these are not measurable for the assembled system.

Gauss method of minimization
The objective function ( 7) is nonlinear, so to find the unknown joint parameters it is necessary to linearize it.Here the Gauss method of minimization is applied [29].First the gradient is found: We search for  = p for which the gradient is zero and denote with ( p) for compactness: However, we cannot easily solve for p, as it appears implicitly in ( p) and ŷAJB ( p).So instead we suppose we have an estimate of p that we denote .We now make two approximations: One, we replace ( p) by (), and two, we make a Taylor series expansion.Using the first two terms of a Taylor series about , we write and the a posteriori term in ( 9) can in the same step be rewritten to: Now (9) becomes This equation is linear in p.For estimates  close enough to the solution, the correction step is small so the linear approximation ( 10) is licit and the new estimate p will be better than  [29].With that, we can now change notation to: p+1 = p and p = , where  + 1 is the number of performed iterations: By rearranging, it is possible to isolate the newest iteration p+1 .By multiplying into the parentheses, taking the terms not depending on p+1 onto the right hand side, and collecting the terms we have: By inverting the left-hand matrix, which we will now denote , onto the right hand side, leads to the final updating expression: where  ∈ R × and  ∈ R  is defined by and  is the number of unknown parameters, +1 is still the number of performed iterations, and ℎ  is a step size control parameter.
Note, ŷAJB  holds the estimated frequency responses for the -iteration found using the joint parameters p .The iterative process lies in the repeated use of (15).At each update the estimate p+1 improves.The minimization procedure is similar to that in [30,31], with the distinct difference that this case takes the adjoint of the matrix instead of the transpose.Fig. 6 illustrates the iterative process for the AM structure with a VP.

Finding the derivative of the frequency-based substructuring expression
To apply (15), it is necessary to determine the derivatives of the responses.The derivative of ŷAJB  w.r.t. the unknown parameters p depends on the case.,internal matches the measured response  AJB internal,mea as best as possible.For better visualization of the process, the illustration is in 2D, but for an actual experiment, the structure is in 3D and the impact points and sensors are placed along all three dimensions.

Case 1
First, a notation clarification for this case.For each frequency line: only given in a vector format instead of a matrix.The derivative of the vector w.r.t.multiple parameters is then a matrix (denoted as  in the previous section) of this form for  unknown parameters.For this reason, it is also so, that the derivative w.r.t.parameter , named , is the same as  ŶAB j   , with the only difference that the first is a vector format, the second a matrix.As the substructuring Eq. ( 4) in Section 2 is in matrix format, it is easiest to determine  ŶAB j   and subsequently reshape it into a vector for the optimization.Differentiating (4) w.r.t each parameter   , for  = 1.. can be written as: as the matrix  A|B is independent of p, this simplifies to: For the sake of compactness, denoting the term to be inverted as , and by following matrix differentiation rules, the derivative of an inverse can be rewritten as: with the differential of  given only by one term since the first term,  A|B  T , is independent of   .This leads to the simplified expression of the derivative of (4): Now, the equation reveals that it is only necessary to determine the differential of the joint's admittance.Note, without the above steps, for systems with a joint of more than two DoFs, it is hardly feasible to find  ŶAB j   numerically, as one needs to determine the derivative of a large inverse matrix (see (22)), a computationally exhausting task.
The last required step is to determine the derivatives of the joint admittance: We know that   is only present in a single entry, Y j  , (as Ŷj is a diagonal matrix per definition -cross-coupling is not possible with the current implementation of the Boolean matrix  and the quasi-static LM formulation).Everything else is zero.So the derivative of the diagonal matrix is: The size of  Ŷj   is  × , where  is the number of joint parameters, but it will always contain only one non-zero entry.For e.g.,  = 1 and  = 6, it is a 6 × 6 matrix with an entrance in (1, 1), and zeros in the rest: We calculate  ŶAB j for each   and reshape subsequently each into the vector . Repeating this procedure for all  will give the matrix  ŷAJB  .For a two-DoF example, thus two interface springs, it is given that: To obtain the derivative of the full system, both are inserted into (25) and ordered as in (20).Note that this expands to more degrees of freedom in the same form.The only limit in this approach is that the joint admittance has to be a diagonal matrix.

Case 2
The derivative of the coupling expression (6) w.r.t. an unknown parameter   is in this case given by: The expression is slightly more complex than the previous case, as several terms depend on   , but by following the same line of thought as before, it can still be simplified to just the differential of the ŶJ itself.First step is to take the derivative of each of the two terms: Then applying the differential product rule for matrices on the second term: Next, is to find the derivative of the last term, again by the product rule: using the rule of the differential of an inverse matrix to find the derivative of the first term above: and differential of  is then given by Going back to the expression (30) and inserting the above relations gives the final expression: In short form: Note, all the differentials have boiled down to  ŶA|J|B , and from Section 2 it is known that only the joint admittance depends on : Here exemplified with a simple two-mass-spring joint (as shown in Fig. 2): As long as there are no cross-coupling terms, the exact expressions above are also valid for the VP interpretations with 12 DoFs, as illustrated in Fig. 3, Case 2. In Case 2, the matrices above expand to 12 × 12 matrices in a block format.E.g., the 2 × 2 stiffness matrix  J is simply repeated along the diagonal, with   ,  = 1..6.Note that cross-coupling can be introduced in the VP, but the stiffness matrix  J will have cross-coupling terms and will be more than just a block-matrix extension of the two-mass-spring joint.In turn, the derivatives of ŶJ will also change, which is why the cross-coupling extension is a future work.
The final step is as for Case 1: The derivatives are inserted back into (37) and ( 35), and finally the vector reorganization is performed.
Note that the procedure easily expands to damping parameters by adding   as an unknown parameter and performing the same expansion.

Algorithm implementation
When implementing the iterative procedure, many decisions must be made.What is the best depends on the system type and noise level.The choices made here ensure the algorithm works robustly for the current experimental setup.Optimization fine-tuning is an endless task, and the choices made here are not necessarily optimal, but they work and demonstrate the strength of the novel method.The subsections below list key considerations and choices.

Considerations when setting the weight functions
The weighting of the measurement data, , versus the weighing of the regularization term  plays a significant part in the final results.Theoretically, when considering the real MAP, the covariances should be used [29], thus  = (cov(  )) −1 and  = (cov( p)) −1 .The true covariance is, however, nearly never available, so when applying the algorithm to a new type of experimental data, it must be expected that some tuning is required to ensure a fair balance between the two terms.A diagonal element in  can be interpreted as the confidence in a measured frequency response at a given frequency line, and  is the confidence in the guess.Normally,  is dominant, as the measurements are trusted much more than the guess of parameters.In this implementation, the weight is also used to normalize the data.The FRFs are normalized via .It is a real diagonal matrix with the same constant entry in all fields: () −1 , where  is the maximum value of the squared response magnitude: and  is a constant weight term, that needs to be tuned to , so there is a balance between the two.Possible future improvements could be to scale the weight function  with the measured coherence of the response.The better the coherence, the more weight is given to the frequency line.Secondly, FRF data is disruptive, changing several factors of magnitude over a frequency range.If the covariance is known, the ratio between  and  is automatically correct in the MAP estimator, and all frequency lines are then scaled and weighted relatively the same.However, if the true covariance is not available,  can be interpreted as a confidence weighting, and then it can be of interest to weight larger amplitudes higher.

The a posteriori estimate -a starting guess
Choosing an appropriate a posteriori estimate is not easy, and the whole point of this work, is to determine these parameters using the iterative process.If it were easy, there would be no reason for implementing any of this.However, the primary purpose of the a posteriori is to regularize the optimization during iterations, thus avoiding negative stiffness and stiffness values of unrealistic orders of magnitude.As a starting guess, one can apply values determined in other works.If the  is not too strongly weighted, the a posteriori term will not control the optimization too heavily.Another option would be to make additional measurements and use, e.g., inverse substructuring results as a starting guess.Inverse substructuring uses the off-diagonal terms of the VPT response matrix.Each off-diagonal entry should, in theory, be straight lines over the frequency range of interest if the stiffness is constant.From [15,22], it is evident that straight lines are not always the case, but inverse substructuring does provide a good idea of the order of magnitude of the stiffness.In the case of damping, a possible starting guess can be to assume proportional damping and use the measured natural frequencies of the assembled system to deduce  and  values, and as initial guess apply that the joint damping is  J =  J +  J , where  J is the initial guess of the joint stiffness.Note, to start the iterative procedure, (15), it is also necessary to set an initialization vector p0 .The initialization plays no substantial role and can be set quite freely.In this work, p0 is chosen randomly from a normal distribution of (1 ± 0.01) p.

Convergence criterion
The iterative procedure needs at least one convergence criterion to stop.The algorithm has converged, when all  parameters   do not change significantly between iterations: where  is the parameter number,  is a small number, and  1 is introduced to avoid the denominator going to zero and is as standard set to 10 −10 .The condition checks if the estimates change between iterations.The criterion must be fulfilled for all .A maximum number of allowed iterations is also set (set to 50 iterations for the numerical validation); if this is reached before convergence, a new a posteriori is chosen, and the algorithm starts over.

Updating the a posteriori during iterations
With the standard MAP method, the a posteriori is chosen beforehand and kept constant throughout the iteration process.It can, however, be useful to update the a posteriori during iterations to speed up convergence and to allow the optimum to be found with gradually less influence of the initial guess.In [30], this is done at every other iteration, and [31] applies a one-time update when errors are below a certain threshold.This work uses a similar approach.When a single parameter is no longer changing significantly, its a posteriori estimate is updated: where  is the parameter number, and   is a relatively small number (though larger than the convergence criterion ).Each a posteriori parameter is updated only once during iterations, and all parameters must be updated before the algorithm checks the convergence criteria.
In some cases, when the measured data is very accurate and the joint model represents the connection to a great extent, it can improve the estimation accuracy further to update the a posteriori one last final time after the convergence criteria (39) is reached for all parameters.So, when (39) is fulfilled for all  (for the first time), the a posteriori vector is updated to p = p+1 .In the numerical validation, this feature is applied.

Validation on the simplest example
An example with a single DoF connection is applied to validate both Case 1 and Case 2. The example is the two 3-DoF systems connected in one point, as shown in Figs. 1 and 2. The unknowns are a spring stiffness and damping coefficient.These will be determined from the measured responses at the internal degrees of freedom and not at the interface.We assume, as would be the real case of an interface, that we cannot measure those particular points.The simulated measured assembled data is  For Case 2, the masses of the joint are set to  1,J ,  2,J = 0.08, and assumed to be known quantities (in order to compare Case 1 and Case 2 fairly).The simulations are performed in MATLAB.The goal is to find values of the connecting stiffness and damping that minimize the difference between the measured response of the assembled system and that obtained by coupling the measured response of the two independent substructures with model parameters.
There are infinitely many options for implementing and deciding on this algorithm.Two different configurations are run to illustrate the effect of these choices.The first configuration a estimates the damping and stiffness in separate optimization steps, and the second b estimates them simultaneously.
The frequency range is 0.01-6 Hz, thus covering the five resonant frequencies of the assembled system.The resolution is 0.01 Hz, giving a total of 1198 frequency lines.The a posteriori values are chosen randomly from a normal distribution between 0.6 and 1.6 times the true value and are not weighted very highly compared to the FRFs.We have that  = (0.05) −1 diag( k−2 J , c−2 J ).The a posteriori guess is constant but updated twice during iterations (see Section 3.3).First, according to (40), and then as described one more final time when the convergence criteria are reached.The convergence and update criteria limits are set to  = 1 × 10 −4 and   = 5 × 10 −4 .The frequency response functions are normalized via .It is a diagonal matrix with the same constant entry .The algorithm will also converge if  and  are proportionally scaled up, and 0.001 should be considered in comparison to the much larger 0.05 in .
In configurations 1a and 2a, when the stiffness is estimated separately, the damping parameter is chosen randomly from a normal distribution between 0.6 and 1.6 times the true value and fixed at that value while the stiffness is estimated.In the next step, when estimating the damping, the stiffness is set to the estimate just found and is kept constant.In principle, this process could be repeated, now using the obtained damping and re-estimating the stiffness, but that is not done here.
To make it realistic, we add complex random noise to all of the simulated measurements  AB J mea ,  A mea and  B mea .Also, without the noise, the algorithm performs perfectly, estimating the true stiffness and damping parameters.The noise is a normal distributed random noise added to the signal, on both the real and the imaginary part, with the levels  real = 5 × 10 −5 and  imag = 1 × 10 −6 (as all other parameters assumed to be in SI units, m/N).Fig. 7 shows an example of one simulation to illustrate the iterative process.The estimated parameters are updated at each iteration, as shown in Fig. 7(d).For each iteration, the objective function  decreases, cf.Fig. 7(c).Figs.7(a) and 7(b) shows one of the 16 FRFs' evolvement over iterations.For each iteration, the response, coupled with the estimated joint parameters, approaches the assembled response  AJB mea , shown with a black line.As described, the a posteriori is updated twice during iterations.The blue and magenta dashed lines indicate the updates at iterations 15 and 21.
To illustrate the enforced noise level in the Monte Carlo simulation, Fig. 8 shows an example set of the  AJB mea frequency response functions.Note that the generated responses of the separate substructures are given the same noise level.
We repeat the process shown in Fig. 7a thousand times for each configuration.Table 2 shows the results of the 1000-sample Monte Carlo simulation.Both the noise and the a posteriori guess change with each new sample.Units are left out, as the example is purely numerical.
First and foremost, Table 2 shows that the algorithm is robust.In this simple case, with this noise level, the algorithm is not particularly sensitive to whether or not the damping and stiffness are estimated together or separately.Despite significant noise, it can accurately estimate stiffness and damping with negligible sensitivity to the a posteriori guess.It works whether we use the weakened formulation (Case 1) or the full coupling formulation (Case 2).Notably, the number of iterations is significantly fewer when estimating stiffness and damping together, making simultaneous estimation more recommendable.The full coupling (Case 2) is more demanding computation-wise (more matrix inversions are needed) than Case 1.With the current implementation, Case 2 takes approximately twice the time of Case 1.
We also test the effect of increasing noise levels.Note that the effect of noise can intertwine with the algorithm settings.As noise levels rise, so does the inherent discrepancy (i.e., regardless of joint values) between the measured assembled response and the coupled response.It can be sensible to relax the algorithm requirements when the dataset is very poor.The algorithm is designed to start over after 50 iterations with a new a posteriori if convergence is not reached.For instance, doubling the noise level (without changing the algorithm requirements) changes the stiffness table values of Config 2b from 100.0 ± 0.32% (Table 2) to 99.85 ± 0.6%.For damping from 1.007 ± 0.91% to 1.02 ± 2.2%.Config 1b results are in the same order of magnitude -still an acceptable level of accuracy for almost all applications.If we track how often the algorithm must restart without converging, that number increases for Config 1b and 2b, but still only a few, ∼15 of 1000, need to restart.In contrast, Config 1a and 2a must restart almost every time, making it much more time-consuming and revealing that estimating stiffness and damping in separate steps is very noise-sensitive and not recommendable.We also test an extreme case and make even the peaks of the responses choppy (tenfold the noise level).We still obtain a reasonable stiffness estimate for Config 1b and 2b (on average ≈2.3% off target with a standard deviation of ≈3%), but the damping estimates fail (on average ≈50% off target).In conclusion, reasonable estimates can be obtained for Config 1b and 2b, even with relatively high noise levels.However, the method becomes less robust, requires longer computation time, and damping estimation fails at lower noise levels than stiffness.

Validation of a 6-DoF connection imitating a VP
We now test the algorithm on a 6-DoF connection to validate the implementation before testing it on real data.It will determine if a more complex system, than the one in Section 4.1, exhibits different estimation abilities for the different configurations.The principle is the same as the previous simple example, but now the unknown joint has six connection points.The simulated system used in configurations 1a and 1b is depicted in Fig. 9.The system is purely theoretical as it is one-dimensional, but as the connection consists of 6 DoFs, the implementation is directly applicable to real 3D structures, where the 6 DoFs instead describe the VP, with its three translational (  ,   ,   ) and three rotational (   ,    ,    ) degrees of freedom.Configuration 2a and 2b are Case 2, where   A schematic of a purely theoretical example for validating the implementation of the 6-DoF joint connection formulation.This case has 6 internal DoFs on both A and B. All internal cross-coupling terms are assumed measured, thus  AJB internal,mea is a 12 × 12 matrix, cf.Fig. 5.At every spring, there is also a viscous damper (left out for compactness).The varying sizes of the mass-boxes do not represent the applied mass values, but are merely to show that they are not the same, and to make it easier to see how the different masses are connected.

Table 3
System parameters for Fig. 9. the joint is a separate substructure with mass and inertia.Table 3 provides the two known subsystem parameters.For Case 2, all the masses of the joint are set to  1−12,J = 0.08 and assumed to be known quantities.The true joint stiffness and damping, used to simulate the ''measured'' responses, are given in Table 4.The simulations are performed in MATLAB.As for the simple example, the responses are contaminated with noise on both the real and the imaginary part, with the levels  real = 5 × 10 −5 and  imag = 1 × 10 −6 .Note that the generated responses of the separate substructures are given the same noise level.Fig. 10 shows a single example of the iterative process for configuration 2b.The estimates jump a little up and down in the first few iterations.The jumping is due to the a posteriori and regulates the optimization so that the estimates stay within a reasonable range; after the first few iterations, the estimates are ''on track'' and converge rapidly towards the true values.

Subsystem
Table 4 shows the overall efficiency of the procedure in the presence of noise.It is evident that despite a relatively high noise floor, the mean estimates are fairly accurate for all configurations.The mean estimates of all the parameters deviate on average less than 2% from the true values.However, the standard deviation in percentage is more than 5% for some parameters of configurations (a).Contrarily, the estimates are almost perfect for configurations (b), with negligible standard deviations.When estimating stiffness and damping separately, configurations (a), it is found that when the stiffness estimates are slightly wrong, then the error is compensated by the following damping estimates; why this approach cannot be recommended unless a further iterative step is introduced (re-estimating the stiffness using the obtained damping estimate).This multi-DoF example clearly confirms and strengthens the result from the simple 1 DoF example that estimating stiffness and damping simultaneously, configurations (b), is a more robust method, and even more so in the case of multiple unknown parameters.Whether it is Case 1 or 2 is still not significant.The methods perform approximately the same.The average number of iterations is also lowest when estimating stiffness and damping simultaneously.It should again be noted that Case 2 is more computationally cumbersome -with the current implementation, it takes approximately double the time of Case 1 to perform one iteration.However, implementing the joint as a separate substructure allows for implementing cross-coupling, which speaks for this formulation.
A general note is that, despite the noise, this simulated system is still ''perfect''.''Perfect'', in the sense that the underlying joint model the algorithm tries to fit the data to is exactly right.There are no cross-coupling springs or anything not taken into account.

Table 4
Mean estimate and standard deviation (std) given in percentage, of joint stiffness and damping for four configurations of the algorithm, with complex noise added to the ''measurements'' (estimates based on a 1000 sample Monte Carlo simulation).Note units are (N/m) | (N s/m).For Config 2 the joint mass is 0.08 per mass (a known value).

Config 1a
Config  For an actual structure, the VP is a simplification, and the joint model is an approximation.Therefore, the accuracy of the regression model is limited by the accuracy of the assumptions.In real examples, it is not obvious where to place sensors and impacts, and there is a risk that the joint is not fully activated in all directions of the 6-DoF VP in the measured data.In the following section, we will try the method on experimental data of a complex structure to see if it has application potential.

Experimental setup and algorithm settings
Fig. 11 depicts the experimental setup in the laboratory.The outer dimensions are 307 × 186 × 40 mm for substructure A and 258×368×20 mm for substructure B. The two components are machined from a single aluminium block.The two parts are assembled with an M10 hex head bolt with a locking nut and top/bottom washers.The structure hung in elastic bands to imitate free-free The separate substructures are also tested, as illustrated in Figs.12(a) and 12(b).The separate substructures  A mea and  B mea are measured near the interface, to be able to make the VP transformation.In the figures, these measurements are the accelerometers and hammer hits (red arrows) close to the connection point.As for the assembled structure, the separate substructures also have sensors and impacts (black arrows) at internal DoFs (same positions).
Regarding the influence of the mass of the sensors: The five sensors used for the internal measurements are placed on the structure of both the reference and the individual substructures; thus, their effect will not influence the joint estimation.The sensors used to perform the VPT are an additional influence.However, the sensors used for the VPT are small, with a weight of 4.2 g each.Altogether, the sensors weigh only 16.8 g, while substructures A and B weigh a hundred times more, 1258 g and 1246 g, respectively.A more problematic issue is the damping coming from their cables.Estimating the damping from the joint is challenging if the separate substructures are more damped than the assembly.Future work could investigate using contactless measurement equipment such as a Scanning Laser Doppler Vibrometer or a High-Speed Camera.For now, we accept that the interface accelerometers introduce an error and add it to the list of bullet points in the Introduction, already counting the VP error, weak nonlinearity, and general measurement uncertainty.
The virtual point transformation of the experimental data is performed using the open-source software PyFBS [32].The iterative estimation is performed in MATLAB.
In the previous section, it was concluded that separate estimation of stiffness and damping was less efficient than estimating these together.Therefore, we test configuration (b) on the experimental data.
For Case 2, we insert a guess of the active mass and inertia of the joint.Looking at the VP substructure in Fig. 3, we need to define the masses  A J, ,  A J, ,  A J, and  B J, ,  B J, ,  B J, , and likewise for the mass moment of inertia .How much mass and inertia to include in this virtual point is not known, but the total mass of the actual bolt and nut (approximately 52 g) can be used as an indicator of the size.After a few tests of the algorithm, we assume the VP connection to be symmetric and that the  and  direction are given by  A J, ,  A J, =  A J, ,  A J, = 15 g.As including the same mass contribution in the  direction leads to divergence of the algorithm (the  direction stiffness converged to a very low value), we assume it to be only one gram at each DoF in .The mass moment of inertia is roughly calculated from the bolt dimensions to 2.2 × 10 −5 N m for  and , and 1 × 10 −6 N m in , around the centre of the bolt.The calculation is made in COMSOL, based on a steel M10 hexagon bolt and nut, by integrating over the domain and the density.
Changing the inertia values to some degree does not influence the procedure, as they are small.In principle, the mass and inertia of the joint could also be treated as unknown parameters and estimated the same way as the stiffness and damping.Estimating also the joint mass is future work.
A few algorithm adjustments were made compared to the simple numerical examples.As the real stiffness is of magnitude 10 4 −10 9 (also varying significantly from rotational to translational stiffness), it is necessary to allow for more flexibility so that parameters can change orders of magnitude; this leads to the risk of negative estimates.Therefore, if one parameter prediction turns negative, the previous parameter value is applied, and the weight of its a posteriori (changing a single diagonal element in ) is doubled to hinder this in the next iteration.Secondly, the final update of the a posteriori, used in the numerical examples, is less meaningful here, as we cannot reach ''perfection''.This step has been removed.
There are many choices to make -should the whole frequency range be used?Should all frequency lines in that range be used?Should some frequency response be taken out due to poor quality?The estimates will change depending on the choice.However, the estimate should stay in the same order of magnitude, no matter the choice.If not, it indicates that this particular parameter is hard to estimate.
The estimates are based on the range 400-2000 Hz, as the first modes are not strongly affected by the flexibility in the joint and since the noise levels were more significant in the first 100 Hz.The measurement resolution is 0.5 Hz, and the estimate is based on all lines.It was checked if it affects the estimate when using fewer frequency lines.If, e.g., only every second frequency line is applied in the optimization.In this case, it does not play a significant role (at least not for the stiffness).An argument for using fewer lines is to speed up the iterative process (as the cost per iteration is less, with less data) to quickly get an idea of the order of magnitude of the parameters.Originally, there was  imp ×  chn = 15 × 17 = 255 measured responses, but for the optimization, 45 were omitted due to bad S/N-ratio.

Results and analysis
Fig. 13 illustrates the iterative procedure for Case 2. As the experimental data is not ''perfect'' and the joint model is an approximation, it requires more iterations to reach convergence (71 in total) than for the numerical examples.The principle is, however, the same.
For each iteration, the joint parameters are updated (cf.Fig. 14) so that the coupling of the experimentally measured FRFs of substructures A and B is more like the measurements of the assembled structure.There are several ways to see that updating the joint parameters improves the correlation.Fig. 13(b) demonstrates it by showing the average coherence (written as COH in all figures), the average Local Amplitude Criterion (LAC) [33] and the average Frequency Response Assurance Criterion (FRAC) [34].The average (red circles) and standard deviations (grey crosses) are based on the 210 measured coefficients to  AJB (45 were omitted due to bad S/N) over the entire range (400-2000 Hz).Using the average coherence and LAC-value as a consistency check is inspired by [35].For reference, the COH and LAC are computed by [35] where * denotes the complex conjugate, and ŶAJB  (  ) is the estimated response at channel  after impact  for frequency   , and Y AJB mea, (  ) is the measured response of the assembled system.The average is used to get an overall consistency.The average coherence is computed as (the same definition for LAC): where  is the number of frequency lines in a response.The Frequency Response Assurance Criterion (FRAC) [34] shares traits with the Modal Assurance Criteria (MAC), made for comparing mode shapes.It directly provides a single number between 0 and 1, representing the correlation over the whole frequency range.FRAC-values are often low [34], as they are sensitive to shifts in the phase.In [34], it is recommended to consider only the magnitude of the response if that is the case.As our interest is especially in capturing the natural frequencies of the assembled response, it is reasonable to make a magnitude-to-magnitude comparison: where  are row vectors describing the response for the whole frequency range at channel  after impact .The average FRAC over all the measured responses is given by: It also appears from , the objective function in Fig. 13(a), that the algorithm works.For each iteration,  decreases.It does not decrease to zero, as the joint model and data are not a perfect match, but it is 32% smaller in the last iteration compared to the initial value.Lastly, the improvements can qualitatively be observed in the FRF plots in Fig. 13(c).Especially at the resonance near 1600 Hz, it is evident that the coupled responses capture that peak better after each iteration.The phase plots are hard to read after 71 iterations, so these are not shown here, but the procedure is, as previously, based on the complex measured responses.
The joint is very stiff, so the standard approach for this structure is to assume a direct VP connection with no flexibility and no joint damping (Case 0).For comparison, such a coupling leads to average COH, LAC, and FRAC values of 0.81, 0.84, and 0.32.The FRAC values should always be interpreted with care, as they are susceptible to the precise position of peaks.However, that is precisely what we want to capture in this investigation.Using a rigid connection (Case 0), several of the resonance peaks are too high up in the frequency spectrum, which means the connection is too stiff.A flexible connection significantly improves this in both Case 1 and 2. The FRAC values are almost doubled to 0.51-0.52.Table 5 shows all the obtained joint parameters and the average correlations for Case 0, 1 and 2.
How to check if the estimates in Table 5 are reasonable?We can check if changing one parameter, a factor of 10, worsens the correlations.If yes, then the order of magnitude of the estimate can be assumed to be reliable.Fig. 15 shows the relative change in average coherence, FRAC, and LAC values when changing one parameter at a time by a factor of 10.Negative values indicate a drop in correlation.How to read the figure: If increasing or decreasing the value of a parameter by a factor of 10 does not influence the correlation (the FRAC, COH, and LAC change is close to 0%), it is not an essential parameter for this structure's global response, and we cannot trust the estimate.However, if the correlations drop when increasing and decreasing a parameter, we can trust that the parameter should be within this range.If there is no drop in correlation when increasing the parameter, we are close to a ''clamped'' parameter.For   that is the case.Here, the estimated stiffness value is so high that it is, in reality, clamped.This also explains the difference between the estimate of Case 1 and 2. One predicts 5.7 × 10 9 N/m, the other 1.3 × 10 9 N/m, but everything in the order 10 9 N/m is effectively rigid.We can also conclude that the translational stiffness in ,  is reliable for several reasons: Direction  and  are both perpendicular to the bolt length, so it makes sense to get approximately the same stiffness in both directions.The same stiffness values are also obtained for Case 1 and 2, and when changing the parameters by a factor of 10, the correlation drops significantly.The rotational stiffness estimates are most likely somewhat unreliable, as there is barely any effect to be seen when changing them.However, it is still useful information that the rotational stiffness in the joint has little influence on the global dynamics compared to the translational stiffness.
The damping estimates are the same order of magnitude for both Case 1 and 2. Notably, both cases in Table 5 give an estimate of  J,  that is two orders of magnitude smaller than the given a posteriori.At the same time, Fig. 15 shows that the average correlation drops when changing  J,  .Both observations indicate it is reliable to assume  J,  is in the order of 1 × 10 1 N m s/rad.The same

Table 5
Experimental estimates of parameters.The noted COH, FRAC and LAC values are the means and standard deviations for all measurements over the frequency range 400-2000 Hz with a resolution of 0.5 Hz.   applies to  J, .The rest of the damping estimates, in  and  direction, do not influence the average correlations very much.To determine how trustworthy these damping estimates are, it is necessary to reconsider how the correlation between measurements and reconstructed dynamics is used in the quality estimation.An average, as performed in (43), can smooth out subtle details.By instead computing the average around a single resonance (e.g., around 1600 Hz), it might be possible to see if changing the damping parameter will influence the correlation at that specific peak.Such investigations will be systematized in future work.

Case
The experimental example demonstrates that considering the joint's flexibility (Cases 1 and 2) leads to better agreement with the assembled reference measurement than ignoring it (Case 0).Both Case 1 and Case 2 perform similarly.Interestingly, including the bolt mass (Case 2) does not significantly alter the estimation, suggesting its negligible influence on this setup.From the numerical validation, we know that Case 2 is computation-wise more costly, at least in the current implementation.However, Case 2 is also easier to adjust to take cross-coupling stiffness terms into account, which makes it the most interesting for further development.As the computation time, for both Case 1 and Case 2, is a matter of a few minutes for this experimental dataset, it is not a weighty argument for choosing Case 1.

Discussion of the method's future application potential
The method presented in this work finds robust and reliable stiffness estimates for a bolted joint.A strength of the method, and also a novelty, is the parametrization before the optimization.The parametrization means that the method directly determines the optimal stiffness and damping rather than the optimal admittance.It can be interpreted as a regularization, ensuring the results stay within a physical solution space.Nonetheless, the parametrization is also simultaneously a restriction.The algorithm tries to find the best parameters but in a limited ''solution space''.In broad terms, it tries to squeeze the problem into a box it does not necessarily fit into.
The current 6-DoF connection model makes the algorithm directly applicable to a range of other structures.Without further adaptation, the algorithm can estimate connection stiffness on any structure with a single-point connection modelled with 6 DoFs (either through a virtual point or directly by connecting six actual measurement points).
Natural next steps can be to expand the underlying joint model to have cross-coupling terms and to add a term in the objective function, ensuring that terms that are not sensitive in the region are forced to be zero to get a minimal number of parameters, despite cross-coupling.Further investigations on how to place sensors on the internal DoFs to ensure the best visibility of the joint parameters could also improve estimations.
At the moment, the weight function of the measurements is size-controlled; each frequency line is weighted w.r.t.its magnitude.In general, trusting large amplitudes more than small amplitudes is not bad, but it is better to weigh the lines around a resonance the highest, even if one peak has a smaller magnitude than another.The goal must be to have the coupling fit the best around the resonances.Secondly, not all resonances activate the joint, so it would also improve estimates to choose to weigh the frequency lines that we expect the most sensitive to the joint parameters higher.However, such steps may be hard to introduce in a generalized manner but are definitely possible for specific structures and examples.
The tests with the experimental structure made it clear that robust damping estimations are challenging to obtain in practice.However, in the current investigation, the evaluation of the estimates' quality was relatively crude.We looked only at average coherences (FRACs and LACs) over the whole spectrum.Damping values will mainly play a role around the resonances in the FRFs, so it would be fair to make a more refined correlation check for future studies.A possibility could be to consider the magnitude around the resonances and ignore the rest of the spectrum.Only the peaks' magnitudes are likely to be affected by small damping changes.
Another work for the future would be to extend the method to include frequency-dependent parameters, for instance to identify joints that exhibit viscoelastic behaviour.The bolted joint of the experimental test structure was in this investigation nominally tightened; but the looser the bolt, the more likely is a non-negligible nonlinear behaviour [26,36].Therefore, it would be interesting to extend the method to nonlinear joints, using, for instance, the Response-controlled stepped-sine Testing method (RCT) approach as proposed in [37] to identify amplitude-dependent properties.
The method's application potential will depend on the ability to make the necessary measurements, and here, the boundary conditions play a role.Our specific experimental setup is free-hanging.In principle, the outer boundary conditions can be anything for the assembled system.However, the two separate substructures must be measured with the same outer boundary conditions.Consider, e.g., a complex example of a car frame and door.Naturally, one will encounter many difficulties when testing such a large and complex setup as a car frame.However, the example can illustrate how to consider boundary conditions.The car frame is substructure A, and the door is substructure B. Imagine that we want to estimate the stiffness of the door hinge.The car door must be measured free-hanging, as its only contact point is the hinge connection.The car frame could have different tyres and be placed on different road surfaces, creating various outer boundary conditions.However, all that does not matter as long as the outer boundary conditions are the same when measuring the separate car frame and the reference assembly (with the door).

Conclusions
The main conclusions are listed below: • A novel method (MAP-FBS) for estimating joint stiffness and damping has been proposed, combining optimization with the physical constraints of an underlying joint model with the knowledge from frequency-based substructuring.• Different versions of this concept have been explored numerically.In Case 1 the joint is quasi-static, consisting only of springs and dampers, and in Case 2, the joint is a separate substructure, including the bolt's mass and mass moment of inertia.Both methods give very reliable estimates when testing numerical data, and they appear equally efficient.• In the current implementation, Case 2 is computationally more time-consuming than Case 1, but estimation time is still only a matter of minutes.Consequently, they are equally recommendable.Which one to use depends on the application (if joint mass is negligible or not).• Different implementations were tested, estimating stiffness and damping simultaneously and in subsequent steps.Numerically, with synthesized noise, simultaneous estimation was the best approach, as it was less likely to diverge and require a restart of the algorithm.When doubling the noise, the estimates are still accurate, but with tenfold the noise (an extreme case) damping estimates fail.• A Virtual Point (VP) connection was introduced to make the method applicable for real bolted structures.Using a VP, the joint connection is simplified to three translational springs ( J, ,  J, ,  J, ) and three rotational springs ( J,  ,  J,  ,  J,  ).Dampers can be added to all six directions.• The novel method was successfully applied to experimental data of a real bolted jointed structure, providing physically reliable estimates of stiffness.The analysis shows the translational stiffness in ,  and  to give the most reliable results.• Further analysis is required to properly evaluate if the experimental damping estimates are reliable since the overall correlation between the responses coupled with the model joint and measured assembled responses is not so sensitive to damping.Nonetheless, the estimates can give an indication of which directions are more damped than others.• The novel method is very general and the current algorithm can be applied directly on any real structure with a 6-DoF connection without further adaptation (e.g. by using a VP to represent a single bolt connection).
Future work will be to test the approach on a range of different types of joints, validate its generality, and test different bolt tension levels to determine the sensitivity of the method.Additionally, the natural next stage is to perform a large comparative study between the different methods listed in Table 1 to determine on which fronts this novel MAP-FBS method is superior/inferior and under what conditions it should be applied.Implementing cross-coupling springs in the Virtual Point could also be interesting, as that may influence some experimental setups.Lastly, developing a more refined way to test the reliability of the damping estimates would be useful.

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.

Fig. 1 .
Fig. 1.Case 1: Conceptual drawing.A linear system consisting of two parts, A and B, connected by an unknown spring  j and damper  j .It should be imagined that the number of interface DoFs can be larger and thus multiple springs can connect the two substructures A and B.

Fig. 2 .
Fig. 2. Case 2: Conceptual drawing.Linear system consisting of three parts, A, J and B. It should be imagined that the number of interface DoFs can be larger and thus a multiple mass-springs system can connect the two substructures A and B.

Fig.
Fig. 3. VP substructure interpretations.Despite the one-dimensional visualization, the VP is a three-dimensional abstraction, providing translational and rotational stiffness in all spatial directions.Note, for every spring there could also be a damper.The superscripts A and B indicates the DoFs for substructure A and B, respectively.

Fig. 5 .
Fig.5.The internal DoFs used for the estimation.The connection for a real structure is always defined as a VP, with 6 × 6 DoFs, but the number of internal DoFs depend on the system.Often is it easier to make many impacts, and only have a limited number of sensors.The number of measurements on A and B is also not necessarily the same.

Fig. 6 .
Fig. 6.The principle of the optimization.The measured quantities are the frequency responses ŶA mea and ŶB mea , at internal and interface DoFs, as well as  AJB internal,mea .By introducing a model ŶJ , and iterating the joint parameters p over  until the coupled response ŶAJB,internal matches the measured response  AJB internal,mea as best as possible.For better visualization of the process, the illustration is in 2D, but for an actual experiment, the structure is in 3D and the impact points and sensors are placed along all three dimensions.

Fig. 7 .
Fig. 7.An example of the iterative process of Config 2b.For each iteration, the coupled responses are updated, illustrated for one response in (a) and (b), the objective function  gets smaller (c), and the estimated parameters are updated (d).In this example the initial a posteriori was kJ = 128 and cJ = 0.66 (∽30−40% different from true values).The convergence criteria is reached after 26 iterations.

Fig. 8 .
Fig. 8. Example of the 16 internal frequency response functions  AJB mea with random noise.

Fig. 9 .
Fig. 9.A schematic of a purely theoretical example for validating the implementation of the 6-DoF joint connection formulation.This case has 6 internal DoFs on both A and B. All internal cross-coupling terms are assumed measured, thus  AJB internal,mea is a 12 × 12 matrix, cf.Fig.5.At every spring, there is also a viscous damper (left out for compactness).The varying sizes of the mass-boxes do not represent the applied mass values, but are merely to show that they are not the same, and to make it easier to see how the different masses are connected.

Fig. 10 .
Fig. 10.Example of the iterative process for the 6-DoF estimation.

Fig. 11 .
Fig. 11.Experimental setup: (a) A and B assembled, (b) substructure A (with zoomed view of the interface) (c) substructure B.

Fig. 12 .
Fig. 12. Overview of measurements, with impact points (arrows) and sensors.In (a) and (b) substructure A and B measured alone.The red sphere in the bolt hole of substructure A marks the position of the VP.Red are the impacts near the interface used for the VP transformation, together with the sensors at the interface.In (c) the assembled structure: measurements are only made at internal DoFs, away from the interface.

Fig. 13 .
Fig. 13.Iterative process for Config 2b, (a) objective function  as a function of iterations.(b) Average FRAC, COH, and LAC values (red circles), and their standard deviations (grey crosses) as a function of iterations, (c) shows the FRFs for a selection of measurements over the iterations.

Fig. 14 .
Fig. 14.How the parameters update over iterations for Case 2.

Table 2
Mean estimate and standard deviation of joint stiffness and damping for four configurations of the algorithm, with a complex noise added to the ''measurements'' (estimates based on a 1000 sample Monte Carlo simulation).For Configuration 2 the joint mass is 0.08 per mass, and it is assumed a known value.Units are (N/m) | (N s/m).