Inverse Problems in Structural EngineeringView this Special Issue
Research Article | Open Access
An Advanced Coupled Genetic Algorithm for Identifying Unknown Moving Loads on Bridge Decks
This study deals with an inverse method to identify moving loads on bridge decks using the finite element method (FEM) and a coupled genetic algorithm (c-GA). We developed the inverse technique using a coupled genetic algorithm that can make global solution searches possible as opposed to classical gradient-based optimization techniques. The technique described in this paper allows us to not only detect the weight of moving vehicles but also find their moving velocities. To demonstrate the feasibility of the method, the algorithm is applied to a bridge deck model with beam elements. In addition, 1D and 3D finite element models are simulated to study the influence of measurement errors and model uncertainty between numerical and real structures. The results demonstrate the excellence of the method from the standpoints of computation efficiency and avoidance of premature convergence.
Recently, moving loads detection technologies have been developed to assess the bridge condition with limited information related to moving loads. In these studies, the measurement of the weight and velocity of the moving vehicle is significant for the design of bridges and pavements, performance assessment, better maintenance, and the control of illegal vehicles on bridges and highways.
A variety of methods have been carried out to identify the weight and velocity of moving vehicles on bridge structures. These methods are called system identification (SI) methods based on the inverse problem because it uses the static or dynamic responses of structures such as natural frequency, mode shape, and time histories of acceleration of the vehicles. At an early stage of the development of the moving load identification method, only static axle loads are measured. Improved methods capable of identifying both the static and dynamic moving loads were developed from the latter half of the 1980s. Hoshiya and Maruyama  identified a moving load on a simply supported beam by applying an extended Kalman filter. O’Connor and Chan [2, 3] developed a method that measured both static and dynamic loads using strain responses obtained from bridge decks modeled as an assembly of lumped masses with the massless elastic beam element. Considering the interaction forces between vehicle and bridge with the viscous damping on the Euler-Bernoulli beam model, Law et al.  identified moving loads using modal superposition principle in time domain. Moreover, Law et al.  proposed a frequency time domain method by performing Fourier transformation of the load-response relationship to identify the moving load directly using least squares. Chan et al.  applied a moving loads identification method based on Euler’s beam theory together with modal analysis. In addition to these methods, Steffen Jr. and Rade  used a Fourier series to identify moving loads on a simply supported beam. To obtain a good quality of moving loads, Chan et al.  developed a method applicable to the prestressed concrete bridges considering the effect of the prestress on the structure. Chan and Ashebo  identified moving loads on a continuous bridge using only a target span. Also, they found that the accuracy of the identified value of moving loads based on bending moment response is better than acceleration response. The effect of the interaction between the bridge and the vehicle, such as the dynamic properties of the bridge and vehicle, made the road surface roughness an important criterion for the identification of the moving load and many studies have been proposed to confirm the results. Cantieni  performed the field test regarding the vehicle-bridge interaction. Hwang and Nowak  contributed to the development of a reliability-based design code by developing a process for the calculation of the dynamic load. Wang et al.  studied the dynamic response of a single span multigirder bridge deck under single and two vehicle loads with different velocities by adopting grillage beam theory and Huang et al.  extended this method to multigirder bridges. Chatterjee et al.  performed a simulation with a quarter truck model by simplifying the continuous bridge deck as a continuous Euler-Bernoulli beam subjected to torsional vibration. In addition, the effect of the speed parameter, the vehicle/bridge frequency ratio, and damping of the bridge and road roughness was reported by Yang et al. .
To solve these inverse problems, direct search methods based on metaheuristics techniques and artificial intelligence, such as genetic algorithms (GAs), simulated annealing (SA) methods, and neural networks (NN), have been introduced and promisingly applied to the field of structural identification . Among them, GAs attract our attention because they do not require a considerable amount of data in advance in dealing with complex problems and make global solution search possible as opposed to classical gradient-based optimization techniques. Suh et al.  presented a coupled neurogenetic technique that is able to identify the region and extent of damage in a beam or frame structure using only frequency information. Mares and Surace  demonstrated the ability of a GA to identify damage in elastic structures. Friswell et al.  combined the genetic and eigen sensitivity algorithms for locating damage. Chou and Ghaboussi  proposed a GA-based method to determine the region and extent of damage in truss structures from the measured static displacements. Krawczuk  presented a wave propagation approach to detect damage in beam structures based on GA and the gradient technique. However, conventional GAs have a limit in solving inverse problems using GA because of the high computational cost of a large number of iterations: it is necessary to perform iterative forward computations for each chromosome. Thus, the total time spent in solving the forward problem could be extremely long, usually on the order of magnitude of several thousand or more depending on the complexity of the problem. Unlike these, a uniform microgenetic algorithm (μGA) can avoid premature convergence and deliver faster convergence to the near optimal region than a simple GA. Carroll  found that a uniform μGA is more robust in handing an order-3 deceptive function than the traditional GA methods. Au et al.  developed a natural-frequency-based μGA for detecting damage in a one-dimensional beam. Lee and Wooh  applied a μGA for detecting damage in plate structures subjected to dynamic loading. However, the numerical examples used in the studies are not practical in that they regard damage as a predetermined rectangular element. Rus et al.  examined the identification of defects in laminated composite structures subjected to in-plane static loads using the boundary element method (BEM). They used an elliptic function with five unknown parameters to represent an arbitrary damage shape. However, this approach has several limitations such as the use of static loads, limitation to in-plane behaviors, and lack of capability to detect the extent of damage. Lee et al.  presented a method to detect stiffness degradations using a modified bivariate Gaussian function, with which it is possible to consider an arbitrary damage shape. The method has been applied to concrete plates subjected to impact loads using the finite element method (FEM) and μGA. Lately, Park et al.  proposed a damage-detection technique to determinate damage distribution, which is a modified form of the bivariate Gaussian distribution function. This method uses a combination of the combined finite element method (FEM) and the advanced μGA and moving loads as input excitation. However, it is required to know the position and velocity of the moving load in advance .
In this study, an advanced coupled GA (c-GA) is developed for solving the inverse problem. The c-GA can avoid premature convergence like μGA and provide faster convergence to the near optimal solution than μGA by introducing a conventional gradient-based technique for which local searching power is sufficient due to fast convergence. In this study, we use a three-dimensional (3D) beam model for obtaining the measurement data. The approach of the present study is more similar to reality than using the same model for generating the simulated experiment and using the inverse procedure. This is because the difference of measured data due to modeling error between real bridges and numerical models can be accounted for. A one-dimensional (1D) beam model is used to run an inverse procedure to obtain better computational efficiency. In this study, the model uncertainty is significantly considered as a realistic circumstance.
2. Forward Procedure
For the transient analysis of a bridge deck subjected to the effects of moving loads, an implicit time integration method, called the Newmark integration technique, is adopted with the integration parameters and , which lead to constant-average acceleration approximation. Considering a moving load with a velocity on a plate element, the total moving distance () of the load at time is given by where denotes the initial coordinate of the moving load in the longitudinal direction. The location number of the element that the moving load passes through at time can be expressed as where and and are the number of division elements in the longitudinal () and transverse direction (), is the initial coordinate of the moving load in the transverse direction, and are the lengths of the plate in the and directions, respectively, and INT() means the integer part of the value in parentheses.
The moving load vectors at an arbitrary location on the th element of the plate should be inevitably distributed into the nodal loads using the zero-order Hermite (Lagrange) interpolation function . The natural coordinates of the element for the moving load at time can be derived as
In a four-node element with three degrees of freedom per node, the moving load distribution into four neighborhood nodes not considering distribution of moment can be expressed as
The total external force vectors applied on the plate at can be obtained by summing the distributed loads as given by
In the Newmark integration scheme, the effective loads at time can be calculated as
The dynamic displacements , accelerometers , and velocities at time can be solved as where the triangularized effective stiffness matrix is and , , , , and are integration constants in the Newmark integration method, respectively.
The governing equation of motion of the system is written in the form where and are the displacement and acceleration vectors, respectively; is the mass matrix without loss before and after damage; is the stiffness reduction matrix; is the time history of the applied moving load.
3. Coupled Genetic Algorithm
The accuracy of a system identification method applying a conventional gradient-based technique may be influenced by the insufficient initial information. That is, the trap of minima due to the incomplete initial values often occurs and this problem should be dealt with. In this study, to resolve the earlier trouble of a gradient-based technique in estimating system parameters, the advanced system identification technique applying the c-GA is employed. In general, GAs are good at global searching but slow at converging because they are necessary to perform iterative forward computations for each chromosome. On the other hand, a conventional gradient-based technique is good at local searching but lacks a global search power; thus, to enhance searching capability and improve convergence performance, the incorporation of a GA with a conventional gradient-based optimization technique is enough to attract attention. The introduced c-GA is such an algorithm. In addition, the c-GA can overcome limitations in solving inverse problems using a conventional GA because the high computational cost of a large number of iterations can be reduced dramatically by operating on a very small population size. The small population size very often leads to the phenomenon of genetic drift in chromosomes over a few generations. To maintain the genetic diversity in the population, the c-GA enhances the genetic operation and search strategy. To identify the unknown parameters, the c-GA should be combined with the FEM that can reflect the change of structural properties and loading condition in bridges in the present state. In other words, the FE model parameters that can explain the change of stiffness due to damage under moving loads sensitively should be utilized as identification parameters in an inverse procedure. In terms of the genetic algorithm, the unknown parameter vector is represented by the selected individuals in each iteration. Meanwhile, dynamic analysis using the combination of FEM and c-GA can be considered the following vector function: where is the vector space of identification variables; is the number of identification variables; is the vector space of dynamic response such as displacements or acceleration data; is the number of measured data used to identify a system; are dynamic responses calculated by using the dynamic analysis combining the FEM with the c-GA from an arbitrary . The function reflects the distribution of stiffness reductions in structures transformed equivalently from dynamic responses changed from damages. Thus, the system identification for computing a distribution of stiffness reduction is described as the following optimization problem set: where are the measured data and () is the function satisfying . Figure 1 illustrates a flow chart for identifying the parameters computed by the combination of FEM and c-GA, as applied in this study. Using the combined finite element analysis and coupled genetic algorithm, the location of a damaged region as well as the distribution of deteriorated stiffness finally can be determined by investigating the unknown parameters .
4. Numerical Examples
4.1. Comparison with Other Algorithms
The performance of the coupled genetic algorithm (c-GA) proposed in this study is compared with that of the well-known simple genetic algorithm (SGA) and the microgenetic algorithm (μGA). The test function used is called the foxhole function and the form is as follows: in which the coefficient is continuous, nonconvex, nonquadratic, and multimodal, has low-dimensional detection property, and, especially, has 25 local solutions. The global solution exists at and the function value at the point is unity. In Figure 2, convergence rates obtained through the SGA, the μGA, and the proposed c-GA are represented. The figures indicate that the convergence rate of the c-GA is faster than the others. The reason is that the c-GA decreases the size of the population by 1/5 compared with the SGA; thus, the computation amount in a forward procedure is reduced, which increases the power for detecting the optimal global solution by introducing the gradient-like selection technique in reproduction operation and adopting enhanced strategies such as the elitist strategy and the scaling windows scheme compared with the μGA.
4.2. Moving Force Detection
4.2.1. Numerical Model
To consider the uncertainties between 2D model and 3D reality, the measurement data obtained from the actual bridge modeled by a three-dimensional FE model [29, 30] shown in Figure 3 are used in the inverse procedure for detecting the characteristics of moving loads. Because of the uncertainties that occurred from the different models, attention is needed to the selection of the measurement range and location for the moving load excitation. In this study, the dynamic responses such as acceleration, velocity, and displacement are measured at the bottom of the 3D solid model and moving loads are excited at the center line on the top of it as depicted in Figure 3. The geometrical and material properties of the beam are the same as those of the numerical model. For the verification of the ability to identify moving loads using the proposed technique combined with the FEM and the c-GA, in this study, three cases of numerical tests are carried out. To focus on the identification of moving loads, we consider an undamaged concrete beam’s 20 divided elements subjected to the moving load of unknown velocity and weight, as shown in Figure 4. The length, height, width, and density of the beams are 24.0 m, 0.8 m, 0.4 m, and 2,400 kg/m3, respectively. In implementation using the coupled genetic algorithm, we examine 4~5 individuals due to moving load excitation cases and the probability of uniform crossover was set as 1.0. The combination of given possibilities can be encoded in binary digits to form each individual. In our case, we use an individual containing 68 chromosomes.
The main goal of this study is to identify the weight and velocity of the moving load passing on a bridge by adopting the quadratic function. In this study, the weight () of the moving load will be considered as a constant and a quadratic function is applied to describe the characteristic on the velocity of moving loads because the actual velocity of the moving load may take shape in a smooth curve and will be changed continuously. Thus, the moving velocity of the vehicle is expressed as where denotes the initial velocity of the th vehicle (or axle); and denote the coefficients of the vehicle with respect to time. Therefore, general identification variables for detecting moving load properties are given by where denotes the number of vehicles considered; denotes the number of the coefficients used to present the characteristic of velocities; denotes the number of the axle weights considered. With different moving load parameters in (13), three examples (EX1~EX3) for four moving load excitation cases (Case I~IV) are given in Table 1.
To quantitatively represent the accuracy of the known parameters identified, the relative percentage errors (RPE) are calculated with respect to the estimated values using the following equation: where and denote the identified and the original values of the known parameters, respectively.
4.2.2. Detection of Parameters
The estimated values of identification variables of several examples (EX1, EX2, and EX3) for Cases I, II, III, and IV based on the c-GA are represented in Table 1 with RPE in parentheses. Figure 5 shows the comparison of the measured velocity and estimated velocity of the vehicle for EX2 in Case I. As the loading time of a moving load increases, the difference in velocities becomes more or less increased but the overall values are found to not be significantly different. The measured and estimated vertical displacements at the midpoint about the same example and case are plotted in Figure 6. The best fitness function values for EX2 in Case I are shown in Figure 7. It is confirmed from this figure that the estimated values of the parameters would be close to real values when the number of generations is increased. The other examples for EX1 and EX2 in Case I showed a similar tendency.
Case II, representing bridges subjected to moving loads of multivehicles, is to describe a more realistic situation about a moving load excitation. Unlike Case I, all moving loads are assumed to move with constant velocities. In addition, two velocities and two weights of the moving loads are estimated. Figure 8 shows a comparison of the measured and estimated velocities of the vehicles for EX1 in Case II. The difference between the measured and estimated velocities is found to be very small. The estimated vertical displacements at the midpoint about the same example and case agree well with the measured ones as depicted in Figure 9. The best fitness function values are plotted in Figure 10. The convergence rate of unknown parameters for Case II is relatively slow compared to Case I since the number of moving loads is double and the moving load properties contained in the dynamic response are superposed and, thus, each property of the moving load is not easily identified. However, it can be observed from the figure that the value of the fitness function converges after approximately the 3,000th generation. The other examples for EX2 and EX3 in Case II showed a similar tendency. In Case III, two axle loads are considered to describe midsize vehicles such as a car. Unlike Cases I and II, in this case, a moving load is assumed to move with constant acceleration. Two axle weights distributed into the front and rear wheels in a vehicle are estimated. Figure 11 shows a comparison of the measured and estimated velocities of the vehicle for EX1 in Case II. The initial velocity of the vehicle shows a slight difference but the overall values are found to not differ significantly. In this case, a vertical velocity is used as the dynamic response data for identifying the unknown parameters. The vertical velocities at the midpoint of the 3D-reality model and the equivalent model are plotted in Figure 12. The behavior between the two models is shown to be identical from the figure. It is seen from Figure 13 that the estimated values of the parameters would be close to real values if the number of generations is increased. The other examples for EX2 and EX3 in Case III showed a similar tendency.
In Case IV, three axle loads are considered to represent large vehicles such as a truck. Like Case III, in this case, the moving load is assumed to move with constant acceleration as well. Figure 14 shows a comparison of the measured velocity and estimated velocity of EX3. As depicted in the figure, the velocity of the vehicle is detected exactly. The quantitative properties of the moving load are also present in Table 1. Initial velocity and constant acceleration are relatively exactly evaluated as 0.224% and 2.120%, respectively, in RPE. The accuracy of estimation results of axle loads is unsteady in this example but, in the remaining examples (EX1 and EX2), for the most, part axle loads are exactly evaluated with the relative percentage errors being less than 5%, as expected. In this case, accelerations are used as the dynamic response data for identifying unknown parameters. The accelerations at the midpoint of the 3D-reality model and the equivalent model are plotted in Figure 15. From the figure, we see that the behavior between the two models is shown to be similar. The best fitness function values are plotted in Figure 16, and it is seen that the estimated values of the parameters would be close to real values if the number of generations is increased.
5. Summary and Conclusions
In this study, unknown moving loads parameters in a beam-type structure are estimated by using the c-GA, which has a superior computational efficiency. The damage characteristics are not considered to investigate only properties of moving loads. A quadratic function is applied to describe the characteristic of the velocity of moving loads. In addition, to consider modeling error in this study, dynamic responses obtained from a three-dimensional FE model under a moving load are used to obtain measurement data. Parametric case studies showed that the proposed technique combining the FEM and the c-GA is adequate to detect the properties of moving loads. Based on the present computational results, the following conclusions may be derived.(1)The c-GA, in comparison with its predecessor (SGA) or other conventional searching techniques, is more attractive not only because it can avoid premature convergence but also because it converges faster.(2)The weight and velocity of the moving load are estimated with small error for examples in Cases I, II, III, and IV considered in this study.(3)The velocity properties of moving loads are relatively evaluated by using a quadratic function for representing continuous moving velocity.(4)The detections based on the dynamic response such as displacement, velocity, and acceleration are all effective in the inverse procedure.
It is concluded from numerical examples that the proposed method works well for the numerical experiments that were tested. However, our moving load characterization using two-dimensional modeling is a limited example for solving the inverse problem of complex structures with large degrees of freedom; hence, more advanced studies should be carried out for more realistic bridge models such as a slab or continuous bridge. In addition, for the application of an uncertain measurement system, it is necessary to develop filtering techniques that are capable of canceling noise signals that may be contained in the measured data. Furthermore, for more complicated situations like damaged bridges subjected to unknown moving loads, more advanced studies that can simultaneously detect the characteristics of the damage and the moving load are needed.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
This research was supported by a basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (no. 2012R1A1A1014722).
- M. Hoshiya and O. Maruyama, “Identification of running load and beam system,” Journal of Engineering Mechanics, vol. 113, no. 6, pp. 813–824, 1987.
- C. O'Connor and T. H. T. Chan, “Dynamic wheel loads from bridge strains,” Journal of Structural Engineering, vol. 114, no. 8, pp. 1703–1723, 1988.
- C. O'Connor and T. H. T. Chan, “Wheel loads from bridge strains: laboratory studies,” Journal of Structural Engineering, vol. 114, no. 8, pp. 1724–1740, 1988.
- S. S. Law, T. H. T. Chan, and Q. H. Zeng, “Moving force identification: a time domain method,” Journal of Sound and Vibration, vol. 201, no. 1, pp. 1–22, 1997.
- S. S. Law, T. H. T. Chan, and Q. H. Zeng, “Moving force identification—a frequency and time domains analysis,” Journal of Dynamic Systems, Measurement and Control, vol. 121, no. 3, pp. 394–401, 1999.
- T. H. T. Chan, S. S. Law, T. H. Yung, and X. R. Yuan, “An interpretive method for moving force identification,” Journal of Sound and Vibration, vol. 219, no. 3, pp. 503–524, 1999.
- V. Steffen Jr. and D. A. Rade, “Identification method of multi-degree-of-freedom systems based on fourier series,” The International Journal of Analytical and Experimental Modal Analysis, vol. 6, no. 4, pp. 271–278, 1991.
- T. H. T. Chan, S. S. Law, and T. H. Yung, “Moving force identification using an existing prestressed concrete bridge,” Engineering Structures, vol. 22, no. 10, pp. 1261–1270, 2000.
- T. H. T. Chan and D. B. Ashebo, “Moving axle load from multi-span continuous bridge: laboratory study,” Transactions of the ASME: Journal of Vibration and Acoustics, vol. 128, no. 4, pp. 521–526, 2006.
- R. Cantieni, “Vehicle/bridge dynamic interaction for highway bridges,” in Proceedings of the 2nd European Conference on Structural Dynamics (Eurodyn '93), vol. 2, pp. 961–968, Trondheim, Norway, 1993.
- E.-S. Hwang and A. S. Nowak, “Simulation of dynamic load for bridges,” Journal of Structural Engineering, vol. 117, no. 5, pp. 1413–1434, 1991.
- T.-L. Wang, D. Z. Huang, and M. Shahawy, “Dynamic response of multigirder bridges,” Journal of Structural Engineering, vol. 118, no. 8, pp. 2222–2238, 1992.
- D. Huang, T.-L. Wang, and M. Shahawy, “Impact analysis of continuous multigirder bridges due to moving vehicles,” Journal of Structural Engineering, vol. 118, no. 12, pp. 3427–3443, 1992.
- P. K. Chatterjee, T. K. Datta, and C. S. Surana, “Vibration of continuous bridges under moving vehicles,” Journal of Sound and Vibration, vol. 169, no. 5, pp. 619–632, 1994.
- Y.-B. Yang, C.-H. Chang, and J.-D. Yau, “An element for analysing vehicle-bridge systems considering vehicle's pitching effect,” International Journal for Numerical Methods in Engineering, vol. 46, no. 7, pp. 1031–1047, 1999.
- T. Marwala, “Damage identification using committee of neural networks,” Journal of Engineering Mechanics, vol. 126, no. 1, pp. 43–50, 2000.
- M.-W. Suh, M.-B. Shim, and M.-Y. Kim, “Crack identification using hybrid neuro-genetic technique,” Journal of Sound and Vibration, vol. 238, no. 4, pp. 617–635, 2000.
- C. Mares and C. Surace, “An application of genetic algorithms to identify damage in elastic structures,” Journal of Sound and Vibration, vol. 195, no. 2, pp. 195–215, 1996.
- M. I. Friswell, J. E. T. Penny, and S. D. Garvey, “A combined genetic and eigensensitivity algorithm for the location of damage in structures,” Computers and Structures, vol. 69, no. 5, pp. 547–556, 1998.
- J.-H. Chou and J. Ghaboussi, “Genetic algorithm in structural damage detection,” Computers and Structures, vol. 79, no. 14, pp. 1335–1353, 2001.
- M. Krawczuk, “Application of spectral beam finite element with a crack and iterative search technique for damage detection,” Finite Elements in Analysis and Design, vol. 38, no. 6, pp. 537–548, 2002.
- D. L. Carroll, “Chemical laser modeling with genetic algorithms,” AIAA Journal, vol. 34, no. 2, pp. 338–346, 1996.
- F. T. K. Au, Y. S. Cheng, L. G. Tham, and Z. Z. Bai, “Structural damage detection based on a micro-genetic algorithm using incomplete and noisy modal test data,” Journal of Sound and Vibration, vol. 259, no. 5, pp. 1081–1094, 2003.
- S.-Y. Lee and S.-C. Wooh, “Waveform-based identification of structural damage using the combined finite element method and microgenetic algorithms,” Journal of Structural Engineering, vol. 131, no. 9, pp. 1464–1472, 2005.
- G. Rus, S.-Y. Le, and R. Gallego, “Defect identification in laminated composite structures by BEM from incomplete static data,” International Journal of Solids and Structures, vol. 42, no. 5-6, pp. 1743–1758, 2005.
- S.-Y. Lee, T. Park, and G. Z. Voyiadjis, “Detection of stiffness reductions in concrete decks with arbitrary damage shapes using incomplete dynamic measurements,” ASCE Journal of Engineering Mechanics, vol. 134, no. 7, pp. 567–577, 2008.
- T. Park, M.-H. Noh, S.-Y. Lee, and G. Z. Voyiadjis, “Identification of a distribution of stiffness reduction in reinforced concrete slab bridges subjected to moving loads,” Journal of Bridge Engineering, vol. 14, no. 5, pp. 355–365, 2009.
- M.-H. Noh and S.-Y. Lee, “A bivariate Gaussian function approach for inverse cracks identification of forced-vibrating bridge decks,” Inverse Problems in Science and Engineering, vol. 21, no. 6, pp. 1047–1073, 2013.
- ABAQUS, “Defining a model in Abaqus,” Online manual, Version 6.9, 2009, http://abaqusdoc.ucalgary.ca/v6.9/books/usb/default.htm?startat=pt01ch01s03aus03.html#usb-int-imodel.
- ABAQUS, Abaqus Analysis User's Manual, Online Manual, Version 6.9, 2009, http://abaqusdoc.ucalgary.ca/v6.9/books/usb/default.htm.
Copyright © 2014 Sang-Youl Lee. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.