Swarm Intelligent Optimization Conjunction with Kriging Model for Bridge Structure Finite Element Model Updating

: For the simple bridge structure, the ﬁnite element model established by drawing and elastic mechanics method is accurate. However, when faced with large and complex long-span bridge structures, there are inevitable differences between the ﬁnite element model and the physical model, where the model has to be updated. It is problematic that the updating structural matrix cannot be fed back into the existing general ﬁnite element calculation software in the traditional structural matrix updating method. In this paper, a parameter-type updating method based on the “Kriging model + swarm intelligence” optimization is proposed. The Kriging model, based on Genetic Algorithm (GA), Bird Mating Optimizer (BMO), and Particle Swarm Optimization algorithm (PSO), is introduced into the ﬁnite element model, updating this to correct the design parameters of the ﬁnite element model. Firstly, a truss structure was used to verify the effectiveness of the proposed optimization method, and then a cable-stayed bridge was taken as an example. Three methods were used to update the ﬁnite element model of the bridge, and the results of the three optimization algorithms were compared and analyzed. The results show that, compared with the other two methods, the GA-based model updating method has the least time due to the small computation. The results of the BMO-based model were time consuming compared to the other two algorithms, and the parameter identiﬁcation results were better than the GA algorithm. The PSO algorithm-based model updating method to solve the ﬁnite element model was repeated, which required a large amount of computation and was more time consuming; however, it had the highest parameter correction accuracy.


Introduction
During the bridge structural analysis, the finite element model of the structure is usually established according to the design drawings. In order to obtain the real dynamic characteristics of the structure, the parameters of the model must be set accurately [1][2][3]. By analyzing the structure with an accurate finite element model, the static and dynamic responses of the structure can be predicted, and the damage condition of the structure can be directly simulated [4,5]. Therefore, it is of considerable significance to obtain the real finite element model of the bridge structure [6][7][8]. However, in general, due to the complexity of the bridge structure itself, it is challenging for researchers to build models that match theory exactly with the experiment.
The errors between the calculated results of the theoretical model and the test results are attributed to three aspects: (1) model structural error, (2) model parameter error, and (3) model order error [9]. The model error is usually related to the selected mathematical model, which depends on the main characteristics of the actual structure. The error of the model order is caused by the finite element discretization, which cannot be avoided in finite element analysis. Therefore, for most of the bridge structure finite element model updating, the key is to use appropriate methods to reduce the second type of error, that is, the model parameter error [10].
The method of frequency response function [20] can provide adequate response information by itself; however, when using this method, response information needs to be obtained in a high-demand environment, which limits the application of this method to some extent. Because of its strong learning ability and nonlinear mapping ability, the neural network is suitable for finite element model updating. However, the method has some defects in robustness. The genetic algorithm [21] is a worldwide optimization search algorithm. It is natural and universal, robust, suitable for parallel processing, efficient, and practical. The selection of fitness function determines the speed and effect of the algorithm. The simulated annealing method [22] is a standard method of probability calculation. Because of its high quality, a strong initial robustness, simplicity, generality, and ease of implementation, it is often used to search for the optimal solution in an ample space. However, this method also has the disadvantages of a high initial temperature, a slow cooling rate, a low termination temperature, and a time-consuming optimization process. The model updating method based on the statistical method only uses the finite element analysis during the initial sample preparation stage. Therefore, it is suitable for all kinds of projects [23]. However, when using this method, feature extraction, parameter screening, and response surface fitting are also needed, and there are problems such as inappropriate feature extraction and multiple types of parameter selection.
Most of the above methods have problems such as too many parameters to be corrected, too much calculation, and low efficiency of model updating. Not all of the parameters to be updated have a significant influence on the structural model. Therefore, the method based on sensitivity analysis should be used to analyze the updated parameters. The problem with sensitivity-based model updating is that the impact of parameters with high sensitivity on the structure may be small [24], while the influence of parameters with low sensitivity on the structure may be critical. Therefore, it is more meaningful to adopt the structural model updating method, which is not based on sensitivity.
Genetic algorithm (GA) [21] is a calculation method based on "natural selection and survival of the fittest" in the theory of evolution. It is a parallel, random, and adaptive search algorithm. Because GA can effectively avoid the problem of local optimization in the search process, it developed as one of the principal artificial intelligence algorithms.
Bird mating optimizer (BMO) is a new kind of swarm intelligence optimization algorithm [25] which has fewer controlled parameters and can avoid local optimization in the optimization process. Therefore, the BMO algorithm has become a lively research direction in structural finite element model modification.
Particle Swarm Optimization (PSO) [26], an evolutionary computing technology developed in 1995, is derived from the simulation of a simplified social model and belongs to a swarm intelligence algorithm.
In this paper, we firstly studied the updating method of the Kriging model [27]. Three kinds of optimization algorithms, GA, BMO, and the PSO, were applied to the updating method of the Kriging model, respectively. Secondly, a truss structure verified the effectiveness of the three methods and compared the time consumption and calculation accuracy. Next, we applied the Kriging model updating method based on three optimization algorithms to a cable-stayed bridge structure for analysis.

Kriging Model
The Kriging model is based on the original data of information samples in a specific region [28,29], and carries out linear unbiased, minimum variance estimation for unknown data with the same characteristics in the selected region. It consists of two parts: the linear regression part and the non-parametric part. The non-parametric part can be regarded as a Gaussian stationary random process, and a polynomial and a random distribution can represent the Kriging model: where y(x) represents the unknown function; β represents the regression coefficient; f (x) represents the polynomial function of the variable x; m is the number of f (x).
In the design space, f (x) can be used to represent the global approximation and z(x) can be used to simulate the local approximation. z(x) is a Gaussian stationary random process, the mean value of which is zero and the covariance non-zero, and it obeys normal distribution N 0,σ 2 . The covariance matrix of z(x) can be expressed as: where R θ, x i , x j represents the spatial correlation function between selected two points from samples, which directly affects the accuracy of simulation; θ represents the parameters of the correlation function; n represents the total number of points in the sample.
To build a completed Kriging model, assuming the test sample points in the region of n × p are x = (x 1 , x 2 , . . . , x n ), where x k i ∈ R p , p is the number of design variables. The corresponding output can be expressed as: β and σ 2 can be estimated as follows: whereβ andσ 2 are estimated value of β and σ 2 , F represents the estimated value vector containing the f (x) for each test point, R is the correlation matrix of the test points, which can be expressed as: To get the values ofβ andσ 2 , the parameter θ should be figured out due to theβ and σ 2 are related to parameter θ. The maximum likelihood estimation is used to solve the minimum value of the following equation: When solving the optimization problem, GA, BMO, and PSO are adopted in the Kriging model to solve the optimization of θ, respectively. After the parameter θ in the correlation function is obtained, the optimal linear unbiased prediction result of the response can be expressed as: where r T (x)α represents the difference value of residual error of the regression function f T (x)β, and the vectorα can be expressed as: r T (x) in Equation (9) can be expressed as the vector of the correlation function between the test point and the unknown point: The Kriging model established can be used to predict the points to be measured. In Equation (3), the simulation of the random process includes a correlation function that can affect the stationary characteristics of the model. Then, express the correlation function of the test point as [30]: where x k i and x k j represent the part k of the test point. Gaussian function can be used as kernel function of correlation function: The steps to update the finite element model by Kriging model are showed in Figure 1 and explained as follows: (1) Determine the variables to be corrected in the finite element model, and then generate a certain number of samples according to the distribution form of the determined variables. (2) Substitute the generated samples into the finite element model of the bridge for modal analysis, and extract the modal correlation information corresponding to the variables. (3) The variables and the corresponding structural mode frequencies are taken as the input and output samples, respectively, and GA, BMO, or PSO is used to optimize the parameters θ in the correlation function of the Kriging model to complete the establishment of the Kriging model. (4) Take the frequency of modal test of the bridge as the input, use the Kriging model to predict the structure modal, and extract the variables by finite element model updating.

Kriging Model Updating Method Based on GA
When dealing with a specific problem, the GA treats the possible solution to the problem as a solution space. The solution space is regarded as a population, and the solution or solution vector is taken as individuals in the population. When the possible solution to a problem is transformed from the solution space to the search space using the GA, the transformation process is regarded as coding. GA randomly generates the initial population and calculates the fitness value of each individual in the population according to the fitness function. Based on the first generation, according to the evolutionary theory of "natural selection, the survival of the fittest", each offspring will inherit a better solution than the previous generation. In the application of the theory of evolution, individuals with high fitness are selected and calculated, utilizing crossover and mutation operators in genetics to generalize new solutions into a new population. Finally, decoding the optimal individuals in the constrained population and the decoded optimal solution can be taken as the approximate optimal solution of the optimization problem [13]. The specific steps of using GA to update the finite element model can be summarized as follows: (1) Code the samples to generate the initial population; (2) Set the fitness function of the population, which is the inverse of the difference between the modal frequency of the model and the measured field frequency. Then the fitness of the individuals in the population is calculated; (3) Conduct individual selection, crossover, mutation to generate new progeny population; (4) Decode the obtained new progeny population and input it into the finite element model to calculate and extract the modal information of the structure; (5) Obtain the frequency error by comparing the finite element modal frequency with the measured frequency. When the error meets the requirement of accuracy, complete the correction process; if the accuracy cannot meet the requirements, then go back to step (1).

Kriging Model Updating Method Based on BMO
The BMO is an algorithm that imitates the mechanism by which birds produce superior genetic offspring. Consider a set of solution in the optimization problem as a bird. In nature, birds are divided into males and females, with females representing the better solution. The mating patterns of female birds are divided into two categories based on the breeding categories of the birds: polyandry and self-mating. The mating patterns of male birds fall into three categories: monogamy, polygamy, and hybridization. In the BMO algorithm, the number of birds in each category is determined by the user, while in the actual situation, the proportion of birds in monogyny and polygyny is relatively large. In contrast, the proportion of birds in polyandry and hybridization is relatively small [25].

Monogamy
If a male x M and a female x i mate, the equation for producing offspring is: where x b represents the offspring; w is the time-varying weighting factor, which decreases linearly from 1.85 to 0.15; r g is the vector of random distribution of elements, and the interval is from 0 to 1. mc f = 0.85 is the variation control factor. r 1 and r 2 represents the random number between 0 and 1. c represents the random number between the number of optimized solutions and 1; u and l represent the upper and lower critical values of the initial solution, respectively.

Polygyny
Polygyny refers to the male birds that mate with more than one female, and the formula for the offspring can be expressed as: where n i represents the number of females attracted by the male birds; x Pg represents the male birds; x ij indicates that the jth female has good genes.

Promiscuity
If the solution of the worst fitness value is mutated then hybrid birds will appear. In the algorithm, a chaotic sequence generates hybrid birds, and a chaotic sequence generates a new feasible solution. The offspring of the hybrid birds are calculated by using Equation (14) with reference to the monogamous birds.

Parthenogenesis
Parthenogenesis birds represent the best result. Calculate the brood as following: where mc f P represents the variation control factor of self-mating the bird population, which increases linearly from 0.10 to 0.75, r 2 and r 3 represents the random number between 0 and 1. µ is the step size and value is 9.1 × 10 −3 .

Polyandry
In the case of polyandry, the female bird will choose several males with good genes for mating in order to obtain offspring with better genes. See Equation (15) for the calculation of the offspring birds.

Steps
In the selection of female or male mates, apply the roulette selection strategy. The higher the fitness value, the greater the probability of being selected. The specific steps of using BMO to update the finite element model can be summarized as follows: (1) Assuming the size of the bird population is m, the probability of selecting the kth bird is calculated as: where f (x i ) represents the fitness function of the problem to be optimized; m is the size of the population.
(2) The selection probability p(x k ) and accumulation probability Q(x i ) of individuals are calculated according to their serial numbers.
(3) A number r = rand(0, 1) is randomly generated between 0 and 1 to determine which interval the number falls in. If Q(x i−1 ) < r < Q(x i ), the interval is selected.
(4) Repeat step 3, and stop when the number of selected offspring reaches the required number.

Kriging Model Updating Method Based on PSO
PSO is a method developed by simulating the predation behavior of birds. Regard the set objective function as the behavior of the birds searching for food. Regard the range of single or multiple variables in the objective function as the flight space of the birds. Regard the space where the solution of the objective function is located as the space where the birds are searching for food in flight. Regard the birds as particle groups and regard each bird as a particle. In the process of searching for and capturing food, carry out information exchange and transmission among individual birds. The flight trajectory and flight speed are adjusted continuously, and then prey is gradually approached. When birds are in the process of finding and capturing food, individual birds continuously share and transmit information. This sharing of information enables each individual to understand their location. The information is then updated continuously to determine whether the target they are locked on to and the flight status is the best combination, and they pass this information on to other individuals. The birds then search near the target and regard the process of finding the target as the process of finding the optimal solution in the optimization problem [26].
In the PSO algorithm, assume there is a group of particles with the total number of m, which belongs to D-dimensional space, and use the position of x i = (x i1 , x i2 , · · · , x iD ), i = 1, 2, · · · m for each particle in this group.
Then set the objective function, namely the fitness function. Through the objective function, the most satisfying position of each particle in the group in the space can be found and then expressed with P i = (p i1 , p i2 , . . . , p iD ). The velocity of the particle is expressed in terms of V i = (v i1 , v i2 , . . . , v iD ). P g = p g1 , p g2 , . . . , p gD represents the optimal position of each particle in the group.
The updated position and velocity of each particle in the particle swarm can be calculated according to the following formula: where x (i+1)d represents the position of the next generation of the particle, which is calculated by the current position of the particle x id and the velocity of the next generation v (i+1)d . w is the weight for speed update, or inertia weight; c 1 and c 2 is the acceleration factor, and they are generally taken as c 1 = c 2 = [0, 4]; r 1 and r 2 represent the random number between 0 and 1. The individual optimal solution for particles can be expressed in terms of P best , and the optimal solution for particle swarm can be expressed in terms of G best . The optimization process of the PSO algorithm can be summarized as follows: (1) It is assumed that the total number of particles in the particle swarm is m, which belongs to D-dimensional space. The inertia weight is w , the acceleration factor is c 1 and c 2 , and the maximum flight speed is v max . (2) Calculate the fitness value f (x) of each particle according to the objective function, compare the current fitness value of each particle with the adaptive value of the optimal individual position, compare the current fitness value of each particle P i with the fitness value of the global optimal position P best , and calculate the velocity v (i+1)d and position x (i+1)d of the particle at the next moment. (3) Iterative calculation. (4) Calculate the fitness value of each particle and calculate the optimal position of the individual. Calculate the optimal position of the population and the optimal fitness value of the population, and update the speed and position of the particle. (5) Determine whether the iteration reaches the maximum number or the convergence threshold. If the iteration reaches the maximum number or the convergence threshold, the optimization will stop. Instead, go back to step 3.

Latin Hypercube Sampling
Based on the theory of small sample learning and prediction, the samples of the Kriging model are limited. The method of selecting sample points is particularly important to carrying out the experimental design. In this section, adopt the Latin hypercubic sampling method for sample experimental design.
Latin hypercubic sampling is a stratified sampling method proposed by Makay et al. [31]. This method can avoid the problem that the sampling points overlap in the local area, thereby ensuring that the sampling points are evenly distributed in the sampling space. Suppose there is a hypercube, the variable dimension of which is n, x i ∈ x i l x i u , i = 1,2 . . . n, x i represents the i th dimension variable. x i u and x i l are the upper and lower bounds of the x i , respectively. The operation of generating m samples using the Latin hypercube sampling method is as follows: (1) Determine the number of samples m to be taken.
(2) Divide the interval of each dimension x i into m non-overlapping cells with the same probability, and divide the original hypercube into small hypercubes with the number of m n . (3) Generation matrix M, whose dimension is m × n. Each column in the matrix M is formed by random arrangement of sequence {1, 2, · · · , m}. (4) Each row in the matrix M represents a small hypercube to be extracted, and a point is randomly selected from each small hypercube to obtain the required m samples.

Verification of Plane Truss
The Kriging model method based on GA, BMO, and PSO algorithms is used to update the finite element model of a 14-span plane truss structure in reference [32], as shown in Figure 2. There are two supports in the model with a fixed support at the left end and hinge support at the right end. Each member in the structure is a round steel tube with an inner diameter of 0.054 m and an outer diameter of 0.085 m. The initial elastic modulus is 210 GPa, mass density is 7800 kg/m 3 , and Poisson's ratio is 0.3. During the finite element model updating, the elastic modulus and density of the material are selected as the parameters to be corrected. Assume that the structural parameters to be corrected obey the normal distribution of X~N (µ, σ). µ is the initial value (i.e., the design value) of the parameters to be corrected. σ = µ × α is the standard deviation of the parameters to be corrected, and α represents the coefficient of variation of the parameters to be corrected. Set the coefficient of variation of the design value of the elastic modulus and density as 5%.
Use the Latin hypercube sampling method to sample the elastic modulus and density, and substitute 100 samples into the finite element model of the bridge structure, and calculate the first six modal frequencies of 100 truss models. Take the first three frequencies of the structure obtained in the model as training samples, and the Kriging model based on the GA algorithm, BMO algorithm, and PSO algorithm to train the samples, respectively, and obtain the corresponding training models of the three methods. According to the measured data, and using the best three groups of a training model to predict, the elastic modulus and density can be obtained. Use the obtained values to input into the finite element model based on three methods and compare with the reference [32]. The specific results are shown in Table 1. As shown in Table 1, the results obtained by the Kriging model based on the GA algorithm, BMO algorithm, and PSO algorithm are close. The maximum modal frequency error of the reference [32] is 2.23%, and the minimum value is −0.11%. The maximum mean error of results based on the GA algorithm, BMO algorithm, and PSO algorithm is 1.69%, 1.28%, and 1.08%, respectively, and the minimum value is 0.00%, 0.12%, and 0.06%, respectively.

Description
The proposed optimization algorithm was applied to the finite element model updating for a cable-stayed bridge. The main bridge is 604 m in length. The half span of the bridge is 13.45 m, which consists of 2 × 0.5 m collision barrier, 2.5 m hard shoulder, 2 × 3.75 m lanes, 0.5 m curb, 0.25 m constant, and 1.7 m half column width. The main bridge is a prestressed concrete cable-stayed bridge with two towers and a single cable plane, as shown in Figure 3.
The cable tower of this bridge adopts the single-column tower, and the height above the bridge deck is 78 m. The cross-section of the tower column is a circular hollow-core. The longitudinal width is 6.9 m, transverse width is 3.4 m, and thickness is 0.8-1.1 m in the transverse direction and 1.5 m in the longitudinal direction. The central piers are round-ended piers with transverse width of 13.8-15 m and longitudinal width of 5.8 m and 7 m, and the thickness is 1-2 m. At the bottom of the pier, a substantial section is used. The foundation is made up of 16 bored piles with a diameter of 2.8 m, and a monolithic reinforced concrete bearing platform with a height of 6 m. The stay cables are dense cable systems and fan space layout with a total of 184. Considering the reference wind speed is up to 61.3 m/s of the design of the main girder, the high strength galvanized steel wire is used to make the parallel steel wire rope. The primary beam cable spacing is 6 m, and the tower cable spacing is 1.6 m. The main beam section adopts a triangular box with a single box and a double chamber. The full width of the box girder of this bridge is 26.9 m, the beam spacing is 6 m, the height of the main girder is 3.2 m, the thickness of the box girder roof is 27 cm, and the thickness of the bottom plate is 25 cm. The cross-section of the main girder is shown in Figure 4.

Modal Test
Modal tests were carried out on the bridge before the traffic officially opened. In order to obtain the vertical modal characteristics of the structure, vertical measurement points were arranged at each cable at both bounds of the main span of the main girder. The whole bridge had 24 measurement points which divide into four groups. Since the bridge was an asymmetrical structure, the sensors were arranged by using the symmetry of the structure when selecting the measurement points, as shown in Figure 5. Figure 5 shows the layout of the 1/4 bridge, which contains two groups of test points. The first group contains V1, V2, V3, V4, V5, and V6 (test reference points), and the second group contains V7, V8, V9, V10, V11, and V6.
By collecting the modal information of the bridge structure, the time-history curves of the structure can be obtained, respectively (see Figure 6). The results of vibration mode and frequency of the cable-stayed bridge were also obtained, as shown in Table 2.   In order to accurately analyze the structure, Midas/Civil software was used to build the model of the cable-stayed bridge, as shown in Figure 7.
In the process of model updating, the selected parameters were elastic modulus E c1 and bulk density γ c1 of concrete of main girder, and elastic modulus E c2 and bulk density γ c2 of concrete of the main tower.
By analyzing the acceleration response of the main girder, obtain the first six orders, as shown in Table 2. The measured modal frequency of the bridge is higher than that of the finite element model, and the error is between 0% and 13%. Hence, the parameters need to be updated to obtain accurate results. This section used the Kriging model to update the finite element model based on the proposed three algorithms.

Results
From Table 2, it is easy to find that the vertical modes of the first order, second order, third order, fifth order, and sixth order of the main girder of the cable-stayed bridge were effectively identified according to the measured data. Therefore, select the frequencies corresponding to the vertical modes of the first, second, third, fifth, and sixth order of the main girder of the cable-stayed bridge as the training samples for calculation and analysis. According to GA, BMO, and PSO algorithms, train the Kriging model after obtaining the optimization results of theta. The measured results were substituted into the trained model to obtain the correction results of the parameters to be updated. Next, evaluate the accuracy of parameter correction by using the 2-norm of analysis frequency error as the evaluation index. The calculated frequency results after the model updating as shown in Table 3. Table 4 shows the modification results of the design parameters to be modified obtained in the three finite element model modification methods.  Table 5 shows the comparison of the time consumed by three methods. By using the same laptop and the same structure, GA methods took half a minute, followed by the BMO method, and finally the PSO method, which took 90 s.
The Kriging model-based model updating method reduced the requirement for sample size in the optimization process and improved the calculation efficiency without sacrificing the calculation accuracy. According to the results of the Kriging model updating method based on the GA algorithm, the frequency error between the training samples in the finite element model and the updating model was the largest, which ranges from 0.15% to 4.39%. When the Kriging model based on the BMO algorithm was adopted, the error range of the results compared with the measured data was smaller than that of the Kriging model based on the GA algorithm, ranging from 0.06% to 4.24%. The result was similar to the updating result based on the GA algorithm. However, in the process of calculation, in order to obtain high-precision results, it was necessary to set a sufficient number of population and number of descendants, which results in the updating method based on the BMO algorithm being more time-consuming than the updating method based on GA algorithm.
Compared with the optimization results based on the GA algorithm and BMO algorithm, the updating results of the Kriging model based on the PSO algorithm have the smallest error range from 0.02% to 3.53% compared with the measured results. Similarly, if high-precision optimization results were needed, sufficient population numbers and the number of descendants should be set in the algorithm, which also has the problem of being time-consuming. Compared with the other two methods, the Kriging model updating method based on the PSO algorithm was the most time consuming of the three methods.
Calculate the Modal Assurance Criterion (MAC) values of the analysis mode and the test mode to illustrate the degree of correlation between the updating model and the test mode. Since the cable-stayed bridge used in this example mainly collected the vertical frequency and mode information of the main girder during dynamic testing, the MAC values calculated here also focused on the relevant information of the main girder, as shown in Table 6. For the PSO algorithm, the MAC is the best, the second one is the BMO algorithm, and GA is the last one.  Table 7 shows the comparison of 2-norm of errors by using different model updating methods. When the finite element model was not updated, the error of 2-norm between the measured frequencies and the corresponding frequencies of the finite element model was the largest. Among the three Kriging model updating methods, the order of the error of 2-norm between the frequency of the updated model and the measured frequency was the GA algorithm, BMO algorithm, and PSO algorithm.

Conclusions
This paper presented a finite element model updating method based on the Kriging model and intelligent group optimization. The method can optimize the parameters of the bridge structure and of updating the finite element model. Based on the proposed optimization algorithm, verify the effectiveness of the model updating methods by using a plane truss supported beam structure. Then, taking a cable-stayed bridge as an example, the proposed method was used to update the design parameters of some components in the bridge structure.
To update finite element model of bridge, using GA, BMO, and PSO methods can obtain consistent results and have good accuracy. This proves that the three methods used in this paper can be applied not only to simple structures but also to complex structures such as cable-stayed bridges.
Among the three updating methods, the PSO algorithm has the highest accuracy, and the GA algorithm has the lowest accuracy. In terms of computing efficiency, the GA algorithm takes the least time. The BMO algorithm requires relatively more sample data in the calculation process. Therefore, the time consumption was increased. The PSO algorithm is the most time consuming because it needs to be solved iteratively. Funding: The first author wishes to thank the support by National Natural Science Foundation of China Grants 51908139, National-local joint Engineering Laboratory of building health monitoring and disaster prevention technology Grants GG19KF003, Department of Transport of Hubei Province Grants 2020-2-5-1, Scientific research project of Wuhan Polytechnic University Grants Y047.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available from the second author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.