Estimation of structural stiffness with the use of Particle Swarm Optimization

The paper presents the theoretical background and four applications examples of the new method for the estimation of support stiffness coefficients of complex structures modelled discretely (e.g. with the use of the Finite Element Model (FEM) method based on the modified Particle Swarm Optimization (PSO) algorithm. In real-life cases, exact values of the supports’ stiffness coefficients may change for various reasons (e.g. order of fastening, state of the contact surfaces, environment changes, etc.). Because of the unknown coefficients, reliable simulations of fixed structure (i.e. mounted, assembled, not in a free-free state) are difficult to perform. The method serves as a tool to obtain good correlation between the FEM of a structure and the experimental data. Simple modal tests are required to estimate the first few modes of the fixed system. The FEM of the structure is considered in a free-free state and the support stiffness coefficients of the FEM are estimated by the proposed method. Further simulations with test-tuned and correlated complete FEM could be performed in time or frequency domain.


INTRODUCTION
The method presented in this paper addresses an important problem of mechanical structures supports' stiffness coefficients uncertainty. Thanks to this method, it is not only possible to achieve better correlation between the results from the experimental modal analysis and structural Finite Element Model (FEM), but also to significantly reduce time required for the elaboration of a model even for complex structures. The support stiffness' coefficients may depend on pretension, fastening order or the quality and state of the contact surface. These details may have a significant influence on the support stiffness' coefficients but including them directly in the FEM may be difficult.
Contact phenomenon is a well-known and commonly investigated problem in modern applications of the FEM based methods. This is obviously a nonlinear problem (Laursen, 2013) and as such, it cannot be directly included in linear dynamic analysis. Of course, linearization of contact analysis results is possible (Choi and Lim, 2004), but variations of the boundary conditions commonly observed in real applications may highly influence the complete structures dynamics. As the result, modes shapes along with natural frequencies may vary (Moon and Shaw, 1983). To avoid the contact problem modal tests are conducted in a free-free state (McConnell, 2006). This allows the FEM to be correlated with the experimentally obtained modal model only in a free-free state, thus neglecting the problem of uncertainty of stiffness coefficients. In most cases, supporting or fixing the structure can be treated as a very simplified problem of adding damping and stiffness elements to the structure (Kaliński, 1997;Wittbrodt et al., 2007). However, even if the support model is simplified, it is necessary to fine-tune its stiffness and damping coefficients to achieve a good correlation (Mao et al., 2010). A model based solely on the geometrical and material parameters of the fixture often leads to unsatisfactory results of dynamic analysis. Selected effects of nonlinear stiffness and damping can be effectively calculated in the frequency domain (Xiao et al., 2017), but this approach is less effective in the time domain. Of course, nonlinear instead of linear techniques of dynamic analysis can be used, but in addition to the obvious high time consumption and high computational costs of nonlinear dynamic analysis, they may also suffer from numerical errors accumulation, especially during long-term simulations. On the other hand, irregularities (e.g. state of contact surfaces, order of pretension) and nonlinearities of stiffness and damping coefficients of the support may have very little or no effect on the quality of long-term simulations, but they may have significant influence on the support stiffness' coefficients, which could be effectively considered as constant in the linearized system. Analysis of the linearized system may be sufficient in many cases (e.g. self-excited chatter vibrations during manufacturing process (Quintana and Ciurana, 2011;Kaliński and Galewski, 2015;Tomków, 1997;Chodnicki et al., 2013).
Direct estimation of the complex and fixed system model can be very time-consuming and in some conditions the results may be unsatisfactory. One of the reasons is high damping values, which strongly influence the reliability of modal tests (Catbas et al., 2004). In most cases, one can quickly estimate only a few modes and such a modal model may be insufficient to investigate specific problems (e.g. problems dependent on computer simulations).
Simplified modal tests can be performed relatively quickly, with limited number of sensors and on the stand with fixed structure. Nevertheless, even a simplified but experimentally obtained modal model can be sufficient to validate a numerical modal model which is suitable for a large number of applications. The proposed method allows for a quick assessment of unknown stiffness coefficients of support, which in turn allows the correlation of the analytical modal model with a simplified, experimentally obtained model. It is based on a linear model in which damping effects are neglected, and only stiffness coefficients are tuned to obtain a well-correlated FEM model.
The proposed method should be considered as a tool to obtain a good correlation of the model, rather as a method of direct estimation of the stiffness coefficients of the fixture (Mao et al., 2010). In the case of over-defined structures with a large number of unknown stiffness coefficients, the insufficient sensitivity of the proposed Particle Swarm Optimization (PSO) cost function and the level of numerical errors may lead to a low accuracy of estimation (e.g. two closely placed supports may be treated as one). However, further simulations may be performed using such an obtained FEM if the support stiffness coefficients are not directly applied. Usually a smaller modal model (which has limited dimensions) than FEM is used, and the supports are only passive elements.
The proposed method of estimation of the stiffness coefficients of the fixture can be considered as a form of the method for updating the finite element model. Unknown stiffness coefficients can be modelled as stiffness coefficients of spring-damping elements (Kaliński, 1997). Model update methods are usually based on experimentally obtained modal model parameters or frequency response functions (Maia and Silva, 1997;Heylen et al., 2007). Updating the finite element model is an important and actual problem (Wang et al., 2018;Zhang and Hou, 2018;Batou and Nabarrete, 2018;van Ophem et al., 2018;Rakshit and Khare, 2019;Chaudhary, 2020;Hu et al., 2021). An analytical approach may require a high quality and complete experimental model which may need to be expanded to allow the least square approach for the FEM update. Insensitivity and ill conditioning are a commonly occurring problem (Kono et al., 2013), even if it is possible to obtain a good experimental model. When analysing the boundary conditions of a fixed structure, Latin American Journal of Solids and Structures, 2021, 18(2), e352 3/18 nonlinearities may also be a problem (Kono et al., 2015). The estimated stiffness coefficients may depend not only on the properties of the physical structures, but also on the values of pretension of the individual supports. Even the order in which those supports have been adjusted can be an important factor (Liu et al., 2019). Moreover, the number of stiffness coefficients estimated may be too high to apply analytical optimization techniques. Due to these problems, the FEM update method was developed based on various forms of Artificial Intelligence (AI) techniques (Maia and Silva, 1997;Levin and Lieven, 1998;Marwala, 2010;Park et al., 2017;Negri et al., 2018;Bartilson et al., 2019;Bianconi et al., 2020;Kang et al., 2020;Naranjo-Pérez et al., 2020;Zhang and Sun, 2020). The main difference of the proposed method is that the computation of the cost function is not based on the modal sensitivity analysis (Fox and Kapoor, 1968;Adhiakri, 1999), but on the low dimensional eigenvalue problem. The reason for this is that the computation of the modal frequency and shapes derivatives may not be accurate for the estimated support stiffness' coefficients which are bounding the structure's FEM, which include rigid body modes. The advantage of the proposed method of estimating the stiffness coefficients of the fixture is the possibility of achieving very good compliance of natural frequencies.
The presented method is based on the author's modification of the well-known PSO algorithm (Eberhart and Kennedy, 1995). Main modifications of the PSO algorithm are the adjusted distribution of the particles to the large scope of the unknown stiffness parameters and introduction of algorithm cost functions to increase performance and avoid local minima. Thanks to the use of the PSO algorithm, very good computational throughput can be obtained because the computation for each individual particle can be easily performed in parallel using independent computing units. In addition, it is now observed that the upcoming generations of CPUs will have a much larger number of processing units, and specialised parallel processors for High Performance Computing are available commercially (Guo et al., 2015). Thanks to this, the application of the presented algorithm will be even more attractive in the near future.

MODAL COUPLING
The theoretical basis of the presented algorithm is the modal coupling technique (Kaliński, 2012). Matrix equation of the dynamics of a linear stationary discrete system in a free-free state, i.e.: can be transformed into modal coordinates: using transformation: where: = [ 1 2 … ] -a matrix of normal modes of a free-free-state system, = [ 1 2 … ] -vector of modal coordinates of a free-free-state system, m -the number of modal coordinates considered.
Modal coordinates in the general case do not have physical sense. By means of an abstract mathematical transformation, they facilitate a significant reduction in the dimension of the dynamic system.
If the matrix of normal modes is M-orthonormalised, then, assuming proportional damping, (2) has the form: where: = [ 1 2 … ] -matrix of angular natural frequencies of the free-free-state system, = [ 1 2 … ] -matrix of dimensionless damping coefficients of the free-free-state system. Equation (4) with additional constraints will be used to compute the angular natural frequencies and normal modes of a constrained system. It is assumed that the number of modal coordinates m contains: -all rigid body modes of the free-free-state system, -a set of the first few elastic modes of the free-free-state system, the number of which should be not lower than the number m c of flexible modes of the constrained system to be computed, -all Guyan modes (Guyan, 1965) of a free-free-state system for all combinations of Degrees of Freedom (DOF) of the supports.
If the proprietary FE software is used to compute the dynamics of a discrete free-free-state system, then a specific option must be used to export rigid body modes and Guyan modes and flexible modes. The stiffness matrix of the linear system with constraints is defined as + T and thus (1) has the form: which leads (similarly to eq. (1)-(4)) to the matrix equation of dynamics of a discrete system with constraints, which is defined as follows (Kaliński, 2012): where: -support stiffness matrix with unknown coefficients where index s denotes that the concerns supports, -coupling matrix, whose values for the constrained DOF (connected to the support) should be 1, and 0 otherwise.
The matrices of masses and damping of supports are neglected here.
The new modal stiffness matrix of the coupled system can be defined as follows: where index ξ denotes that the is the modal stiffness matrix of a complete system with dimension [m×m]. The eigenvalues λ ξ of the matrix are squares of the angular frequency ω c of the complete system with constrains. The c index denotes that the system is complete with constrains. The corresponding eigenvectors T ξ can be multiplied by the matrix of normal modes of the free-free-state system, to compute a first few normal modes of the constrained system: The computed angular frequencies ω c can be directly compared to the angular frequencies ω x measured during the modal tests, where the index x concerns the reference values obtained from the experiments. To compare i th computed vector with the j th experimental vector of normal modes, Modal Assurance Criterion (MAC) can be used (Allemang, 2003): MAC values close to 1 for the normal modes indicate a linear dependence of mode shapes. Even if the MAC factor is generally considered not to be a perfect tool for measuring the correlation between mode shapes (Allemang, 2003), it is used here for fast correlation and is assumed sufficient.

STANDARD PARTICLE SWARM OPTIMIZATION ALGORITHM
Particle Swarm Optimization is an AI evolutionary computational technique, suitable for many different classes of problems (Banks et al., 2007;Banks et al., 2008), including the estimation of modal parameters and the update of FEM (Galewski, 2016;El-Kafafy et al., 2015;Shabbir and Omenzetter, 2012). It was first proposed by Eberhart and Kennedy (Eberhart and Kennedy, 1995) and was inspired by the natural behaviour of flocking birds. The basic form of the PSO algorithm is easy to implement and uses simple rules. However, it often effectively outperforms some other AI algorithms (Elbeltagi et al., 2005;Beheshti and Shamsuddin, 2013). The PSO algorithm solves a problem by moving a number of possible solutions called particles in the search space. Each particle stores information about the current and best solution found so far, and also has the knowledge about the best solution found by all particles from the whole swarm. The movement of each particle is described by the velocity and position update equations that consider the abovementioned information. The steps of the PSO algorithm are as follows.
1. Arbitrary selection of parameters ϖ , φ p and φ g which appear later in the velocity update equation, where index p denotes that the parameter is used in particle's local context, and g -in a global context. 2. For each particle i : • initialization of the position x i (where index i denotes the particle number) by means of a uniform distribution of a random vector of size d within the boundary of problems; • initialization of particle velocities v i for each element of the x i vector; • remembering the current position of the particle as the best known particle position p i = x i to date.
3. Finding the particle with the best value of the cost function f(x i ) and remembering its position as the best known global position g = x i . 4. For each particle i , repeat until the stop criterion is assured: • for each dimension d of x i generate random numbers r p,d and r g,d (∈<0;1>, uniform distribution); • update particle's velocity: • update particle's position: • if f(x i ) has better value than f(p i ) (lower or higher, depending on the problem formulation as it can be minimization or maximization) update the best particle position p i = x i ; • if f(x i ) has better value than f(g) update the best global position g = x i . 5. After stopping the algorithm, g represents the best solution found.
Parameters , φ p and φ g are arbitrarily selected, usually based on experience or simulations, and influence the magnitude of the velocity update step size and the proportion of how much the particle takes into consideration the local and global best results during this update. In the case of problems with a significant number of local optima, the efficiency of the algorithm can be reduced. In general, the main problems in the performance of the PSO algorithm are slow convergence and premature convergence (finding the local optimum). The most common techniques of solving these problems are: adding additional elements to the update equations or repositioning the particles when the algorithm slows down or approaches the local optimum (Banks et al., 2007;Banks et al., 2008;Aote et al., 2013).

THE PROPOSED ALGORITHM
The direct application of the PSO algorithm to estimate stiffness coefficients of supports may not produce good results for the following reasons. - With six DOF on the FE support element, the number of all unknown stiffness coefficients x i may be too large to search the entire space with unknown coefficients simultaneously.
-The range of the stiffness coefficients can be very large.

-
The computation time required to solve the eigenvalue problem of the stiffness matrix K ξ of the coupled system (7) usually dominates the total time required for calculations made for each particle. So, the number of particles should be limited.
-Local minima and large spans can lead to undesired, very high particle velocities.
-A cost function based solely on frequency compliance and mode correlations can lead to an erroneous order of modes and the achievement of one of the local minima.
To solve these problems, the new algorithmic cost function was proposed, and the PSO algorithm was modified accordingly. In addition, several procedures are proposed to address the local minima problem.

Description of the algorithm
The proposed algorithm is based on the standard PSO algorithm presented in the previous section with some original modifications. Modification to the standard PSO algorithm will be discussed at the beginning, and later the algorithmic cost function f(x i ) will be presented.
It is assumed that limits (boundaries) for each x i stiffness (unknown) coefficients are defined by the user. Thus, the algorithm must check whether the future position of the particle (11) is within the given limits. If not, the particle "bounces" in the opposite direction (by changing the sign of the particles velocity).
Due to the large number of supports and the wide stiffness range for each DOF of the supports, it is usually not possible to process all unknown stiffness coefficients x i simultaneously. Thus, the initial stiffness value must be chosen for each unknown. The user can enter the initial values as the best-known initial global position g.
For the selected number of unknown N u coefficients, local iterations are computed (steps 1-5 of the standard PSO algorithm). Global iteration is the implementation of all local iterations of the algorithm for all possible combinations of the sets of unknown stiffness coefficients. The minimum number of local iterations is arbitrarily chosen.
Only a few N u unknowns are processed while computing the local iteration. A suboptimal solution can be found for a selected set of unknown stiffness coefficients of the finite elements after computing this step. As a part of each global iteration the presented algorithm should achieve improvement, minimizing the value of the cost function for the global coefficients f(g), and thus ultimately heading towards the global minimum. However, if the improvement of the result (i.e. cost function value) for the global coefficients f(g) is not satisfactory, then the number of unknowns processed during each local iteration may be increased or another procedure for the local minimum problem may be engaged. Several proposals for such procedures will be discussed later in this work.
The limited number of particles must be adjusted to the (usually very large) span of unknown coefficients. The commonly used initial random distribution of the particles position x i may not be the best solution. Instead, the particles are initially distributed uniformly in a linear manner for a range of low stiffness (e.g. from 1 N/mm to 1000 N/mm). For a larger range of stiffness values (e.g., more than 1000 N/mm), the particles are initially distributed in logarithmic manner. In addition, the best global position g is attributed to the initial particle distribution.
A large span of unknown coefficients and the presence of local minima can lead to too high velocities of particles v i (10). It has been observed that particles can achieve velocities high enough to "jump over" the global optimum solution during a single local iteration. To limit this negative effect, the maximum absolute velocity limit has been introduced for each computed particle velocity v i,d (10) as a factor ε of the actual particle position x i,d : The appropriate factor ε is simply a global value given arbitrarily (e.g. ε = 0.05).
The cost function f(x i ) is based on the compliance of the selected number of reference modes N ref . The first modes of the structure usually consist mainly of Rigid Body Modes. They are more sensitive to changes of support stiffness coefficients values. On the contrary, higher order modes are purely flexible modes, so they are much less sensitive to changes in support stiffness values. It is usually very easy to get a very good correlation with the higher order modes and at the same time get stuck in the local minimum. To avoid this problem, computations start with a minimum set of reference modes (usually only the first identified mode is processed in the beginning) and additional modes are added if the result of the cost function for the global coefficients f(g) is lower than the arbitrarily given threshold. In addition, based on the results of the cost function for the global coefficients f(g) and the number of processed modes, different parameters of the proposed algorithm can be tuned, e.g. weighing coefficients of the cost function, expected frequencies compliance or MAC threshold values.
The partial cost function for the selected reference mode consists of two parts. The first one describes the correlation of mode frequencies, and the second -the correlation of modes: where: ω c , ω x -computed and given (reference) angular frequency, , -computed and given vector of normal modes, 1 2 , w w -weighing coefficients, whose values can be tuned according to the frequency range, units and etc.
Index j denotes that the value of f j is calculated for the j th experimental mode. It must be noted, that considering only the frequency compliance is not sufficient. It is easy to get very good accuracy in the compliance of natural frequencies and a completely uncorrelated set of modes. The correct selection of the computed and reference mode for comparison (13) is the key for effective applications of the proposed algorithm. When tuning the stiffness coefficients, the order of the modes may change (Gallina et al., 2011). Therefore, it seems that it is more important to track part of the MAC value than to achieve first of all frequency compliance, although care should be taken and appropriate weighing factors (w 1 , w 2 ) can be adjusted when computing the partial cost function. The computation of the cost function f(x i ) is based on the selected number of reference eigenfrequencies N ref and is carried out as follows.
1. The reference mode is selected, and the high penalty value is set as a result of the partial cost function.
2. Checking if all obtained eigenvalues are real numbers.
3. The control of compliance of the eigenfrequencies is performed: the difference between the computed eigenfrequencies and reference ones should be smaller than the desired threshold value (reference eigenfrequency multiplied by the precision parameter, which is automatically adjusted during the operation of the algorithm; first, its value is high for the initial state and for example 10.0 can be given, but then it is gradually reduced to 0.1. (9) is calculated for all computed eigenfrequencies.

MAC factor
5. The result of the cost function for the selected reference mode (13) is computed for the best MAC value (close to 1.0).
6. Steps 1-5 are repeated for all actually processed reference modes. The value of the cost function f(x i ) is the sum of the cost functions f j (x i ) for the selected j th reference modes.
In conclusion, the PSO algorithm has been modified by: -introducing the maximum particle velocity, -changing the initial distribution of random particles to linear-logarithmic distribution in order to cover a large span of unknown coefficients, -introducing the algorithmic, nonlinear cost function f(x i ), with variable global parameters defining the precision requirements for frequency compliance and the number N ref of the reference modes considered.
Both numbers (frequency precision and the number of modes considered) are associated with f(g) and change online, during the operation of the algorithm.
In addition, suggested "rescue" procedures (i.e. to escape from local minima) have been defined to solve the problem of local minima. Initialisation procedures have also been defined to quickly evaluate the initial coefficients. The problem of reaching the local minimum is detected when the cost function for the best coefficients f(g) tends to approach high values, and its improvement observed during at least two global iterations is smaller than the arbitrary set threshold. The scheme of the whole algorithm is presented in Fig. 1. Note that some of its details, which have not been described so far, are discussed in the subsequent chapters of this paper.

Local minima avoidance
To detect a situation, where the algorithm is stuck at the local minimum, the actual J and previous J old values of the cost function for the best coefficients f(g) are compared after each global iteration. If the value of the cost function for the best coefficients J min = f(g) is not small enough, and after the next global iteration the observed improvement is not promising (i.e. J−J old > J min ), one of the following rescue procedures can be used. However, if the value of the cost function for the best coefficients f(g) is less than the specified threshold value J ref , the number of reference modes should be increased.
The number of the processed unknowns at each local iteration can be greater than one. Processing two unknown values at once usually solves not only the majority of local minima problems but is also good for quick evaluation of initial values for all unknown coefficients. The disadvantage of this method is the high computing cost. The number of particles and the number of local iterations should be selected to cover the range of unknown coefficients being processed. However, the number of particles increases geometrically with the number of unknown coefficients processed simultaneously. The total number of particles should be adapted to the number of the CPUs cores, because all computations for each particle can be performed in parallel at each local iteration. For example, in the case of two unknown coefficients processed simultaneously 8 x 8 = 64 particles can be used for a quad-core CPU. A larger number of particles requires significantly more computational resources and usually takes more time to perform each local iteration. In addition, the number of local iterations required to complete one global iteration increases as the total number of unknown combinations of coefficients increases. However, it was observed that quite good results can be achieved if a relatively small number of iterations were used (limit was set to 4). In the case of processing more than one unknown coefficient at once, the number of iterations was automatically increased if the results of the cost function f(x i ) were smaller than the arbitrarily set threshold value at each subsequent local iteration. For a specific combination of unknown coefficients, it was worth spending much more time searching for a better solution. The maximum number of different unknown coefficients processed simultaneously was never greater than three. Moreover, if only the improvement of the cost function for the best coefficients f(g) was better than the arbitrarily set threshold J min after each global iteration, the algorithm switched back to searching for one unknown coefficient at a time. Specific thresholds can be set individually for each case considered. If computational resources are limiting factors or the algorithm coincides with the local minimum, a different rescue solution is proposed. During the FEM processing, it was observed that the sums of all estimated stiffness coefficients of all supports for each particular direction were usually well consistent with the sums of all the stiffness coefficients defined in the reference FEM for each specific direction. However, some of the estimated factors have far too high or too small values. Thus, one of the proposed rescue procedures is to use average values for each specific direction as the new initial (best known) coefficients g.
Finally, it is best to avoid local minima problem at all. This can be achieved by providing good estimated limit values for stiffness coefficients or good initial (best known) values of coefficients g. Therefore, it is recommended to first analyse the physical properties of the specific support. For example, it is not wise to adopt high values of stiffness coefficients for a rotary DOF for a spherical support or low stiffness for a translational, massive, solid support.

SUMMARY AND PRACTICAL RECOMMENDATIONS
Experimental modal tests may be quite basic, however, to achieve the best results, the following recommendations are advised, i.e.: for complex shaped structures, vibrations should be measured along each direction of translation in at least several points not laying on the same straight line; the most important are the first modes from the point of view of the FEM correlation; usually the more modes are identified, the better accuracy will be achieved; the normal modes should be well decoupled, so it is usually better to consider more DOF.
The number of reference modes should be limited. Higher order modes are purely flexible and are not as sensitive to changes in the support stiffness coefficients values as the first few modes.
The total number of particles and the number of local iterations should be adapted to the scope of all processed unknowns during the local iteration. If the number of particles is large, the total number of local iterations can be small and vice versa. It is possible to track changes in cost functions f(x i ) to adjust the number of local iterations during computations.
Each local iteration requires solving the eigenvalue problem for each particle. This is usually the most computationally expensive part of the proposed algorithm. However, all computations for each particle during each local iteration are completely independent and thus they can be performed in parallel. It is recommended to adjust the number of particles to the number of processing units used for the computation. For a practically unlimited range of unknown coefficients (e.g. for linear stiffness from 1 to 1x10 12 N/mm), the implementation of the proposed algorithm on a quad-core CPU was the most effective with a small number of particles (e.g. 12) and relatively large number of local iterations used (e.g. 400).
The most important factor affecting the good performance of the proposed algorithm is the use of an efficient algorithm for fast solving of eigenvalue problem (Cosnuau, 2014). Parallel computations were successfully implemented using the OpenMP programming interface. Several procedures (e.g. matrix multiplications and MAC) have been optimized to achieve good performance. In the case of a complex system and poorly adjusted initial stiffness coefficients a long computational time may be needed. For example, with a 15-support system, it took about 24 hours to find the right solution using a quad-core Intel Core i7-6700 CPU. On the other hand, for an aluminium plate mounted in a vice, computational time was a fraction of a second at the same CPU. A direct comparison of the execution times of the algorithm is not possible due to the random nature of the particle's motion and the possibility of a better or worse choice of initial coefficients. For example, if the tuning of stiffness coefficients for one part with unknown coefficients took a few hours, then it took just a few minutes to tune the coefficients for different fastenings of the same part if the previously obtained results were used as starting coefficients (even if the natural frequencies were different for both fixings for about 10 Hz). Thus, for example, it is possible to efficiently apply the proposed algorithm for a series of similarly but not identically fixes parts, to independently correlate the model for each case. It's worth noting that this algorithm should scale very well in computer systems with more CPUs. Multi-core processors and coprocessors are currently available (Jeffers and Reinders, 2013;Intel, 2019a;Intel, 2019b), so an acceleration of the algorithm by using more particles and fewer iterations may be expected.

ILLUSTRATIVE EXAMPLES
This part presents examples of the application of the proposed algorithm. First of all, two different FEMs are used to verify the accuracy of the proposed method. Then the results of two experimental applications were discussed: the first one is a plate mounted in a vice and the other one is a complex structure with a multi-point support.

FEM simulations -a simple case
The FEM model of the plate presented in Fig. 2 was created for the purpose of simulation. FEM of the structure consists of 8-nodes cuboid finite elements. The support is modelled as one 6-DOF spring-damping element (Kaliński, 1997) that is connected to the structure by Multi Point Constraints (MPC) (white rectangle in Fig. 2) to minimize the influence of local deflections in the FEM. The stiffness values of the support are defined in Table 1. The plate size was 175×50×5 mm and it was modelled as made of S235JR steel: Young moduls E = 210 GPa, Poisson ratio ν = 0.3, density ρ =7.95 g/cm 3 . The total mass of the plate was 348 g. Additionally, masses of two accelerometers (20 g each) that were later (see section 6.3) attached to the plate during real modal tests were included in the FEM model as MASS elements to better imitate real dynamic conditions of the considered object.
The number of outputs (placement of sensors) has been limited to 2 (Fig. 2 -black rectangles), and all translational operations (in the X, Y and Z directions) have been processed for each output. In addition, the normal modes of the model in free-free state have additional data for the DOF of the support (MPC driving point). The M-orthonormalised modes were determined using the Intes PERMAS solver (Kaliński, 2012) for both free-free and constrained (reference) models. A total of 22 modes (including Rigid Body Modes and Guyan modes) were computed for the model in free-free state and as a reference -10 flexible modes were computed for the constrained model. For the purpose of the simulation test of proposed method, only 6 flexible modes of the constrained structure were used as a reference, i.e. 1 st , 2 nd , 3 rd , 5 th , 6 th and 7 th . The weighting coefficients w 1 , w 2 (13) were set to 100 and 1000, respectively, arbitrary threshold J min (Fig. 1) was set to 1x10 -8 , and arbitrary threshold J ref and the number of processed reference modes N ref were initially equal 1 and when the number of reference modes was increased, the arbitrary threshold J dim was automatically reduced to 0.5. PSO parameters (10) , φ p and φ g were set to 0.25, 0.1 and 0.3, respectively. The coefficient ε determining the maximum particle velocity (12) was set at 0.05. The initial values for all DOF of the support (i.e. spring-damping element) have been set to 10 4 , the maximum values for the translational DOF have been set to 10 8 and for the rotary DOF -10 10 ; minimum values have been set to 10 for all DOF.
For a selected set of reference modes, the obtained compliance between the reference model (given values of stiffness coefficients) and the estimated, by the proposed PSO-based algorithm (computed stiffness values) was very good ( Table 2). The calculated MAC values were close to one, and the frequency differences were significantly small for the selected reference modes. The proposed algorithm converged within a short time of 20 s. The computed values of stiffness coefficients (Table 1) were similar to the reference values.

FEM simulations -a complex case
The FEM model of the part presented in Fig. 3 was created to simulate the CNC face milling process. The FEM of the structure consists of 249 078 10-nodes tetrahedral finite elements. The original part was restricted by 12 supports, but for the purposes of the computer simulations, the support system has been simplified and limited to only 4 supports. Each support is modelled as a linear 6-DOF spring-damping element (Kaliński,1997), which is connected to the structure by MPCs to minimize the influence of local deflections in the FEM model. For each support, the same values of reference stiffness coefficients were defined (Table 3). Spheroidal cast iron EN-GSJ-400-15 was used in model as the workpiece material (in compliance with the original part): Young moduls E = 169 GPa, Poisson ratio ν = 0.275, density ρ =7.3 g/cm 3 . The external size of the workpiece is aproximately 1050×580×255 mm. The total mass of the workpiece was 160 kg.
The number of outputs (placement of sensors) has been limited to 16 (Fig. 4) and signals towards all translational DOF (along axes X, Y and Z) have been processed for each sensor. In addition, the normal modes of the model in the freefree state have additional data for each DOF of the support (MPC driving point). M-orthonormalised modes were determined using the FEM Intes PERMAS solver (Kaliński, 2012) for free-free model. A total of 54 modes (including Rigid Body Modes and Guyan modes) were computed for the model in free-free state and as a reference -30 flexible modes were computed for constrained model. All reference modes were used during the simulation of the proposed algorithm. Both weighting coefficients w 1 , w 2 (13) were set to 1 and other parameters were set as in the previous simulation example (section "FEM simulations -a simple case").
For all of the reference modes, a very good compliance was obtained (Table 4). The calculated MAC values appeared close to one, which evidences good modal assurance between the analysed systems. The computed frequencies have become almost the same compared to the appropriate reference values. The computed values of stiffness coefficients (Table 3) were similar to the reference values as well. However, due to the greater complexity of the model, the proposed algorithm needed as much as 24 hours to achieve the required convergence.

Simple experimental case
In this example, measurements and modelling were carried out for simplified modal tests of a steel plate mounted in a heavy vice (Fig. 5). The exact dimensions and material data are in this case insufficient to build a good quality FEM model (Fig. 2) mainly due to the uncertainty in the contact area.
As can be seen in Fig. 6, the estimated Frequency Response Function (FRF) and computed stability diagram using the poly-reference Least Squares Complex Frequency-Domain (pLSCF-D) method (Guillaume et al., 2003;Kaliński, 2012) can be easily interpreted, so the experimental modal model is very reliable. The vice was modelled as a 6-DOF spring-damping element (Kaliński, 1997), and its two stiffness coefficients were estimated using the proposed algorithm in a fraction of a second. The other four stiffness coefficients were assumed to avoid singularity because measurements were made only in one direction for which the first two bending modes were computed from the experimental data. A very good correlation was obtained (Table 5). It should be emphasized that the purpose of the considerations presented in this example was not only to identify the values of support stiffness coefficients, but above all to obtain a compatible calculation model, which was confirmed by comparing the computed results (i.e. normal modes, natural frequencies) with their reference values.

A complex experimental case
The workpiece model (shown in Fig. 3) has been modified to include all twelve supports (Fig. 7). All supports were modelled as 6-DOF spring-damping elements (Kaliński, 1997) in such a way that there were totally 72 unknown stiffness coefficients. Modal tests of the real structure (Fig. 8) were carried out for a limited number of measuring points (Fig. 7). In this case, 16 piezoelectric accelerometers were used to measure the acceleration in directions perpendicular to the surface in each accelerometer mounting point. Using the pLSCF-D method (Guillaume et al., 2003;Kaliński, 2012), only 7 modes can be identified (Tab. 6). They were used as reference modes for the proposed method. It was assumed that only poles with pole lines reaching the order of the model below 20 in the stability diagram ( Fig. 9) are considered stable. Also some modes must be ignored due to the high non-collinearity of their phase angles. Due to the complexity of nonlinear supports, they can be observed during experimental investigations, and therefore the correct identification of all modal parameters may be limited. Similarly to the simulations, flexible, rigid and Guyan modes were computed using the Intes PERMAS solver (Kaliński, 2012), although different MPC definitions were used in this case (due to the greater number of supports). Correlation results are presented in Tab. 6. As one can see, not all computed modes were identified during the experiment. For some of them (e.g. for the frequency of 430.6 Hz), the computation modes do not find their counterparts in the reference model due to the very low MAC (0.01). Nevertheless, the correlation results are good enough for the FEM model to be later used successfully during computer simulations of the milling process.
Such simulations were performed using a model obtained through the indirect computation ( (6) and (7)). The reason for this was not only that the indirect calculation model did not differ much from the model directly computed using proprietary FE software, but mainly due to the computation time, which can be significantly shorter. In addition, proprietary software was used to compute models with different stiffness coefficients, and simulations were carried out at the stand in an industrial environment. In case of simulation of the milling process, accurate matching (error less than 1 Hz) of natural frequencies is usually more important than perfect MAC values (0.85 and above should be sufficient), because the accuracy of the natural frequencies can influence not only the quantitative results of the simulations, but also its quality. For example, development of self-excited chatter vibrations during machining may be related to the correlation between spindle speed and the natural frequency of the object being manufactured. Thus, proper frequency identification is crucial for process parameters choice (Kaliński et al., 2019;Kaliński et al., 2020).

SUMMARY AND CONCLUSIONS
The presented results of computer calculations and experimental studies allow to formulate conclusions regarding the effectiveness of the presented method of estimating the support stiffness using the modified PSO optimization in order to obtain the expected correlation between the computed modal model and the reference modal model obtained from experiments. The proposed method is an effective solution to the problem of fine-tuning the support stiffness coefficients. The proposed method based on PSO optimization with some original modifications resulted in improved accuracy and, above all, reduced time needed to obtain the expected solutions. This is very important because achieving good accuracy of natural frequencies only using FEM can be very difficult and time consuming. The operation of modal models with dimensions much smaller than the corresponding structural models of objects is also not without importance due to the significant reduction in the demand for memory resources required during the computation process.
The proposed approach allows achieving very accurate compatibility of natural frequencies. This may be of key importance for selected applications (e.g. tool-workpiece vibration analysis, where the effectiveness of testing the stability of the machining process may depend on the precise identification of the frequency values).
After several tests of the program on experimental examples, the authors began to use it almost routinely not only for the purpose for which it was developed, but also in many other cases. Since this helped to obtain a good correlation of the FEM test relatively quickly, it was noticed that it is a very useful tool in scientific research and industrial applications.