Abstract
Load identification, or input identification as the more general term, is a field of study that requires a wide set of disciplines, which suffers from uncertainties caused by the challenges within each discipline. When making load identification, several different approaches exist. For all (or at least most) methods, however, some sort of system model is required. This model may be simple or complex, depending on the system at hand. Typically, if the identification process is vibration fed, the system model will be created from modal parameters. These parameters, however, are often subject to uncertainty and thus may be considered as stochastic variables. In this paper, the root causes of uncertainty for load identification are demonstrated using classical identification techniques. From a numerical perspective, uncertainty is quantified through Monte Carlo simulations. Two results are outlined: one where the identification process is completely blindfolded in its most naive form, and one where the spatial distribution of the load is predefined. In general, it is found that fixing the spatial distribution of the load can compensate for truncation errors in the modal parameters.
1. Introduction
Modal analysis is a convenient and efficient method for retrieving the dynamic properties of a civil engineering structure in a condensed form. Traditionally, two methods exist for modal analysis—experimental modal analysis (EMA) and operational modal analysis (OMA). The return from a modal analysis consists of the following modal parameters: mode shapes, natural frequencies, and modal damping [1]. Modal analysis is a commonly used approach, yet the confidence in the extracted modal parameters will depend on the quality of excitation and response measurement. Although new methods are continuously surfacing to overcome some of the challenges, a remaining uncertainty on the estimates persists. In the application for civil engineering structures, another range of challenges may arise, for instance, when the stationary assumption is violated through variations in environmental conditions such as temperature or mass loading.
The exertion of indirect load identification may be a successive step to the modal analysis. If the retrieved modal parameters accurately and in whole can describe the dynamics of the system, they can be used in the inverse calculation needed when doing indirect load measurements. The principle is that if a dynamic model is available (from EMA or OMA) and the response of the system is measured, then an estimate of the load can be obtained through deconvolution. Load identification has been in focus in the recent decade, and many attempts have been made in order to obtain a stable and accurate method. In recent years, identification techniques in real time using Kalman filters have proven successful in certain cases [2–7]. Besides the recursive model description, another benefit from these methods is the capability to merge different sensor type information. For these filters to perform well, however, several noise models must be determined, which may be difficult. Aucejo et al. [8] thoroughly examined how the different fitting parameters may influence the performance of identification though they omitted model errors from their analysis. Using identification techniques in the time domain, Wang et al. [9–11] studied the interval envelope on the load estimation given uncertainty in the response measurements and integration errors.
Another approach is to deconvolute the response in the frequency domain [12–15]. This was the main approach from the early days in the automotive and aviation industry [16, 17]. The method is based on establishing the relation between the system response and force at different frequencies. The technique requires a recording time, which is why the method does not perform well in real-time applications. Comparative studies between the two approaches are also seen in [18]. The response measurements are often recorded in terms of displacements, velocities, accelerations, or strains. These can be obtained through attached sensors (invasive) or from noninvasive techniques such as lasers, digital image correlation, or even by acoustics [19, 20].
For every method referred to above, a system model is needed. This model may be derived from a finite element (FE) model or from modal parameters directly obtained from the EMA or OMA. Now, an obvious question arises: how accurate must the modal parameters be in order for the load identification to be satisfying? Although some authors already have addressed some of these questions [21–24], this paper will revisit some of the fundamental challenges from a visual point of view. By studying a practical example, the sensitivity of estimation errors will be demonstrated through a Monte Carlo simulation. In order not to disappoint the reader, we note that this paper deals merely with the consequence and not mitigation of the estimation errors in modal parameters.
2. Simulation Setup
A plane cantilever beam will host the basis for this study. The beam is made from Euler–Bernoulli beam elements with two degrees of freedom (DOFs) at each node—a translational and a rotational as indicated in Figure 1. We assume that the beam is proportionally damped through the Rayleigh coefficients α and β. The mass is assumed to be distributed along the elements. The beam is discretized by three elements and fixed at the bottom node which yields, in total, six free DOFs and hence six modes of vibration. We note that the FE model has not yet fully converged for all modes at this nodal resolution; however, this is not of importance for the given study.
The beam model is then subjected to a loading history from which the dynamic response is simulated at a sampling frequency of 1 kHz. The loading time history is selected so that both dynamic response and quasistatic response will be visible within the output. The load is composed of a square impulse followed by a low-frequency wave, again followed by a superimposed high-frequency wave train (15 Hz). Figure 1 shows the system in focus for this current study. The four-noded cantilever beam constitutes the system model. The loading history (input) is shown in the left-hand side of the figure. Each line corresponds to a load time series for the corresponding node, i.e., a flat line resembles zero loading. The solid lines are the horizontal force, whereas the dashed lines are the moment. The loading is set to act at one node only—the second node from the top. The right-hand side shows the corresponding response (output) for each node. Solid lines are the horizontal displacements, and the dashed lines show the rotation of the nodes. Amplitude scaling is irrelevant as the system is linear.
The process of load identification will hence be the reverse order of the process shown in Figure 1. First, the response of the structure is measured. Then, if the dynamic properties of the system are known, we can make an inverse calculation and thus estimate the load causing the measured response. Before proceeding to the load identification, a small amount of noise is added to the system response and hereafter considered as the measured signal (100 dB Gaussian white noise is added).
3. Load Identification Process
In order to assess the sensitivity of the load estimate given the variations in the estimated modal parameters, a method must be chosen. For this paper, two methods will be outlined—both operating in the frequency domain. We have selected two methods that easily can be applied using the results from an EMA or OMA (Typically, if modal parameters are obtained from an OMA, problems persist as the modal mass and hence the mode shapescaling is unknown. For now, we will ignore this challenge and assume that the mode shapes obtained are correctly scaled). The methods originate from the automotive and aviator industry [22, 23] and are often referred to as the “transfer path analysis matrix inversion” or just “pseudoinverse technique.”
3.1. Method 1
The first method, referred to as Method 1, is the most simple and straight forward. This method might be a tempting first approach, but it shall be demonstrated later on why it could lead to unfortunate results. Consider the response of a linear, time-invariant system which, in the frequency domain, is given as follows:where is the system response, is the frequency response function (FRF), and is the load, all given at the discrete frequency ω. This process resembles very well what has been illustrated in Figure 1. The FRF matrix can be computed as a sum of N modal contributions as follows:where is the mode shape vector, is the mode scaling function, and is the complex pole, all for mode r. is the transpose, is the complex conjugate, and is the Hermitian transpose. See [25] for the derivation of this formula.
Now, we presume that the output is available through measurements—both the translation and the rotation of all DOFs. All that is left to do is to calculate the inverse of the FRF and premultiply on both sides in order to obtain an estimate on load .
Since errors or modal truncation can cause the inversion of the FRF to become singular, the inversion is performed using the Moore–Penrose pseudoinverse [26]. Method 1 has been used in [15, 27, 28] with the minor change that the inversion of the FRF matrix is performed through a singular value decomposition technique. Regularization techniques have also shown to reduce the sensitivity in the matrix inversion [29].
3.2. Method 2
Another method, Method 2, will be added for comparison. This method resembles Method 1 to a high degree with only one minor, yet important, difference. In Method 1, there were no restrictions to the solution, and the load estimate could be distributed to any DOF in the system. If the spatial distribution of the load is known (and unchanging), this can be incorporated into the identification process through a separation of variables as follows:where constitutes the spatial distribution of the load and is its scaling function. If the given loading scenario allows this separation and if is known, equation (1) and (3) can be rewritten aswhich then can be solved for the scaling constant :and the final load estimate then reads
Methods to derive the distribution in a similar system have been shown in [13, 30–32], but for this study let us assume that the distribution is a known quantity; hence, .
The equations, (1)–(7), describing Method 1 and Method 2 directly yield a frequency domain estimate of the load. All the following figures in this paper will show the load estimate as a time history in line with the input (left-hand side) in Figure 1. The equivalent time-domain representation is obtained by means of the inverse Fourier transformation.
3.3. Condition of the Frequency Response Function
The condition number may be a measure of how sensitive the matrix is to inversion, where a large condition number indicates an ill-conditioned matrix and vice versa [33]. Singular value decomposition may separate the FRF matrix into singular vectors and singular values. Letting the condition number be defined as the ratio between the largest and smallest singular value, we thus obtain a quantitative measure on the sensitivity. Since the FRF matrix, is defined for a range of frequencies, each frequency will be inherent to a condition number. The condition number(s) for the full rank FRF matrix is shown in Figure 2 as a function of frequency. Here, it is seen how the condition number increases at the resonance frequencies.
There is no sharp limit to assess whether a condition number will yield a stable inversion, but a condition number larger than 100 should call for immediate concern which is the case for most frequencies in our case [33]. Since the FRF matrix is constructed from a set of modal contributions, the rank of the FRF matrix will depend on the number of independent modes. If modal truncation is present, this causes the FRF matrix to be rank deficient, which means that the smallest singular value becomes close to zero and consequently yields a condition number going to infinity.where m and n are the number of rows and columns in while and indicates the condition number and rank, respectively.
If the FRF matrix is constrained by a spatial distribution vector as it is done in Method 2, the product reduces to the vector . Singular values of an nonnull vector is a single nonzero element arranged in an vector [34]. This consequently provides a condition number of unity, which means that the “inversion” will be stable for all frequencies.
4. Modal Truncation
The number of DOFs in a real structure will, of course, be approaching infinity, but the modes which can be identified are depending on the sensor distribution and capability. Consequently, truncation of the modal space will always be present in real-life situations, which is why this section is included [35, 36].
From equation (2), it is seen that the response of a system is a linear contribution from all modes in the system. If modes are omitted from this sum, this will truncate the response function. Figures 3 and 4 show how this modal truncation affects the estimate on the load estimate from the inverse calculations. One additional mode is removed for every case. All modal parameters are treated as deterministic and at their true values, i.e., the only error in the FRF is the higher modes being omitted.
As seen in Figure 3, the load estimate is acceptable when all six modes are used despite the small content of noise on the response measurements. However, when omitting modes from the FRF matrix, the load estimate from Method 1 quickly becomes erroneous, both in terms of distribution and magnitude. Even though the resonance frequency of the modes omitted may be far from the frequency of the load, this truncation still makes a great impact. Meanwhile, in Figure 4, it is seen that Method 2 retrieves more consistent estimates on the load. Any DOF, besides DOF 3, obtains a zero-load estimate following the definition of the spatial distribution given in equation (7). For both methods, it is seen that the estimate for the square impulse overshoots at the discontinuities caused by the Gibbs phenomena for the truncated modal space [37].
If we consider the resynthesized response from the load estimates from either Method 1 or Method 2 and the truncated response function , we observe
Here, ε is an arbitrary error function between the measured response and the resynthesized response . The solution from Method 1 will yield a better response approximation than Method 2. This means that Method 1 is a more mathematically accurate solution. Yet, when comparing Figures 3 and 4, it is clear that in a physical sense Method 2 is more consistent. This is one of the major challenges when making load identification solely based on the response using a least-square approach. Consequently, the load estimate may return as an equivalent loading rather than the actual, as seen in Figure 3.
For the given simulation, the load is conveniently acting at a point where the response is being “measured.” In other cases one might not be this lucky and expansion is needed in order to estimate the response at locations that were not recorded originally. Several expansion techniques (/virtual sensing techniques) exist for this [38, 39], but effectively, similar truncation errors will be introduced as no new modes are added during the expansion.
5. Stochastic Modelling
The modal parameters needed to establish the response function may be obtained through an experimental identification process—either EMA or OMA. The parameters obtained in this process will always be subject to variations. The variations can originate from physical causes such as influence from temperature, operational mass loading, or other nonanticipated nonlinearities violated by the model description. Also, nonphysical uncertainties related to data postprocessing and pole extraction are known to exist [40–43]. These magnify if the excitation of the structure is unfortunate or if the sensor resolution and location is poorly chosen. Noise and other limitations on the sensor may also be a cause of uncertainty.
The variations in the modal parameters can be implemented from many different uncertainty models. The most common is probabilistic, where a probability density function is used to describe the statistical variations. Other formulations such as the fuzzy-set model or a nonprobabilistic interval method could also have been utilized. For this study, however, we will rely on the probabilistic approach and assume that all these aforementioned causes of uncertainty effectively can be modelled as an uncorrelated stochastic variation in each of the estimated modal parameters. The uncertainty in natural frequencies is assumed to follow a normal distribution given a coefficient of variation; hence, estimates for the lower frequencies are the most certain. The coefficient of variance is chosen as for all frequency estimates. The stochastic model for the natural frequency of mode r hence readswhere the expectation μ is assumed to be in line with the true natural frequency for mode r, i.e., . The standard deviation is derived from the coefficient of variation as .
Damping estimates are known to be the most uncertain parameters to quantify. Since the damping typically involves fairly low values—say a few percents—and the estimates are subjected to a vast amount of uncertainty, it must have some skewness in order to avoid negative estimates. The study [44] also noted a positive skewness in damping estimates, yet the distribution fit was not studied. For this simulation, we assume that the distribution can be modelled as a standard gamma distribution with shape parameter α. The stochastic damping parameter model for mode r hence readswhere the shape parameter is assumed to be equal to the true damping ratio for mode r (in percent). That is, which yields a mean .
We assume that mode shapes always will be real following the proportionally damped system. The variations in mode shapes are implemented as a stochastic process described by a spatially uncorrelated variation in individual DOFs and a variation in scaling of the mode shape. For mode r, this meanswhere is the true mode shape for mode r, is a diagonal matrix whose entries contain the uncorrelated noise on individual DOFs, and κ is the scaling of the mode shape. Both noise models are assumed to be normally distributed with a mean μ and standard deviation σ as
The mode shapes are assumed to be mass normalized. Thus, the alteration in mode shapes will affect the modal mass as well. The spatially uncorrelated alterations may violate the mode shape orthogonality, but only to a small degree as the error standard deviation for the mode shapes’ DOFs is set at 0.5%. The consequences of omitting modes have already been shown; hence, only small variations are now included.
Using the principles of Monte Carlo simulation, a realization of the stochastic modelling is shown in Figure 5 for the first three and the last modes. The uncertainty band for the mode shapes is shown using cubic interpolation in order to show the rotational DOFs. The modal assurance criterion (MAC) [45] between the true mode shape and the stochastically altered is shown for each mode. We note that the variations in MAC values are small compared to what may be experienced from experimental work.
6. Results
Before turning to the identification process, let us have a look at the FRF matrix. The stochastic variables shown in Figure 5 are fed into equation (2) through a Monte Carlo simulation, and a sample of the FRF matrix is shown in Figure 6. For the 100 000 Monte Carlo simulations, the upper and lower 98% quantiles are indicated in the figure as the hatched green area. The black solid line is the mean value. The consequences of missing modes in the FRF matrix were studied by Maes et al. [46] and shall not be repeated here.
In order to evaluate the sensitivity in load estimates, the frequency response function is now considered as a stochastic process given the natural variation in modal parameters. Keep in mind that the “measured” system response is kept the same for all simulations and at its true value with a minimum of noise added. In the following three sections, the parameter variations are introduced a little at a time so that it will be more clear what happens to the estimates. Since the mean value of the modal parameters is in line with the true values, the averaged value of the estimates converges towards the true value as well. All six modes will be used in the making of the FRF matrix.
6.1. Natural Frequencies
First, only the natural frequencies are treated as stochastic variables, and the remaining modal parameters are kept as deterministic and exact. For the 100 000 Monte Carlo simulations, the FRF is synthesized, and a load estimate for each simulation is obtained through both methods. Figure 7 shows the upper and lower 98% quantiles for the time-domain load estimates. It is seen that both methods yield a reasonably stable result despite variation in natural frequencies. Three points in time (a, b, and c) are extracted to highlight the error distribution on the estimate. The relative error for the two methods is given by the error standard deviation and summarized in Table 1.
6.2. Natural Frequencies and Damping Ratios
Next, in addition to variations in natural frequency, the damping ratios are also treated as stochastic variables. That means only the mode shapes are left deterministic. The variation in the estimated load is shown for the two methods in Figure 8. As expected, the uncertainty bound for the estimates increases as the damping ratios are also treated as stochastic parameters. Method 1 seems to outperform Method 2, given this level of uncertainty in the modal parameters and noise in the output signal. A ringing effect in the estimate is observed at the discontinuities near the square impulse. In Method 1, however, it is seen how the load estimates are spreading to different nodes to compensate for the errors in the FRF matrix.
6.3. Natural Frequencies, Damping Ratios, and Mode Shapes
Finally, all modal parameters are now considered as stochastic processes, and the consequent result for the estimated load is given in Figure 9. Although the mode shapes are altered by values down to one or two percent for each node, the uncertainty in the load estimates for Method 1 shows an exponential growth. When mode shape errors occur, Method 1 fails to predict the load distribution and consequently the scaling of the load. For reference, the sample FRF shown in Figure 6 hosts the basis for the load estimate given in Figure 9. Note that in Figure 9, the error probability density is unequally scaled for Method 1 and Method 2.
Each of the load estimates shown by Figures 7 to 9 has three highlighted points in time, where an error probability density is shown. The corresponding standard deviation (std) is given in Table 1. We note that Method 1 seems to perform better than Method 2 as long as the mode shapes are intact. However, when introducing minor changes to the mode shapes, Method 1 cannot fully predict to which node the load is applied, which results in a poor estimate on DOF 3.
7. Discussion
We have studied the performance of the two methods given the natural variation in the modal parameters used for the model description. The load estimates when using Method 1, shown in Figures 3 and 9, are not as diverse as they might appear at first glance. If the resulting force is considered instead of every single entry in the estimated load vector, this leads to a more appealing result. First, reexamining the estimates shown in Figure 9 and by summarizing the contributions from each DOF, an equivalent baseline load is obtained. This resulting load estimate is compared with the true baseline load and shown in Figure 10. Now, it is seen how the uncertainty band is narrowed down. For this static system, an equivalent global loading can be estimated using Method 1 despite the stochastic variations in modal parameters as long as all modes are represented.
(a)
(b)
Returning to the modal truncation example in Section 4, we once again compute the resulting load estimate from Method 1 as previously shown in Figure 3. Using the same analogy as explained above, a resulting baseline load is estimated for different number of modes included in the FRF. As the estimates from the truncation study are deterministic, we choose to plot each result on top of each other in Figure 11 together with the true baseline load. We see that the overturning moment is fairly accurate regardless of the number of modes included, while the shear force is drifting when modes are being omitted. When modes are omitted, Method 1 fails to estimate the point of attack and thus the scaling of the load, as it is seen in Figure 11. It should be noted that the degree of overshooting at the impulse is increased as modes are being omitted. Whether the resulting moment estimate remains stable is likely dependent on the static system.
(a)
(b)
There are some items which have not been covered by this study and which should be mentioned as they limit the conclusions.(i)We have chosen a simple static system with plane deformation. The system is blessed by not having any closely spaced modes and is behaving perfectly linear. Now, if the system is more complex, which may be the case for most civil engineering structures, the generalization of the number of modes needed and the precision on the modal parameters may differ.(ii)We have not touched upon multiple load sources or a moving load.(iii)For this paper, the response function is driven by displacement measurements. It has not been assessed whether acceleration or strain measurements would have yielded different results.(iv)Only one noise level on the output signal has been studied.(v)From Figure 8, it may seem like the variation in damping ratios and frequency yields nearly no effect on the load estimate. If the system is subjected to a harmonic load at a frequency near a natural frequency, these modal values may be vital for the load estimate.
8. Conclusion
The sensitivity in load identification following the uncertainty in modal parameters has been studied for two different methods. The uncertainties have been introduced stepwise into the system to demonstrate the influence of different parameters. It was found that fixing the spatial distribution of the load can compensate for the modal truncation in the response function. If one mode is dominating in the response, parameter estimation following this mode is crucial for the following load identification. However, modes with a resonant frequency above the frequency of the load still contribute to the quasistatic response, and omitting these modes will cause an error. A possible idea is to compensate through static deflection shapes/Ritz vector in the FRF matrix, but this has not been included. The spatial distribution of the load was fixed in Method 2 through the distribution vector . If the assumption regarding the load distribution is not correct, this will introduce some systematic errors in the load estimate. Demonstrating this has been omitted from this paper.
Whether or not the stochastic modelling outlined in Section 5 is appropriate for a real-life structure will be left for the reader to decide. However, it has been demonstrated how variations in these parameters may alter the estimates on the loading, and shown how a sensitivity analysis may reveal flaws in the algorithm.
As a general observation, when doing vibration-based load identification, the algorithm for estimating the input must be supported by additional information besides the system response. This additional information may be in terms of the spatial distribution of the load—which was shown with Method 2—or it may be other load models which are driven by wind speed measurements, wave gauge readings, local pressure measurements, or possibly information about the frequency content of the loading. Without aiding the inverse calculations with one of these, the uncertainty of the estimate will exceed what is acceptable. For most of the successful methods available in the literature, the spatial distribution of the load is also defined as a fixed measure, which leaves room for further research in this field.
Appendix
A. Examining Different Loading Positions
In order to demonstrate that the sensitivity is not uniquely related to the position of the load, a few additional cases are included in this appendix.
A.1. Case I
For the first case, the load is moved to the top node at DOF 5 (see Figure 1). Any other settings are the same as described previously. We jump to the result, where all modal parameters are treated as stochastic variables. The results are shown in Figure 12, and again three selected points in time are highlighted. The error standard deviations are given in Table 2. Note that for Method 2, we once again assume that the distribution is correctly foreseen i.e., .
A.2. Case II
Next, we examine the consequences, if the loading is applied as a moment instead. Here, the moment is applied in DOF 4, i.e., at the same node as the base case shown in Figure 1. Again, we jump to the stage where all modal parameters are treated as stochastic parameters. The results are shown in Figure 13, and the error standard deviation for points a, b, and c is given in Table 3. For Method 2, again we assume that the distribution is correctly foreseen, i.e., .
B. Evaluating the Number of Sensors
The final case is dedicated to the study of how many sensors are required for a successful estimate. It has already been shown that for the truncated modal space, Method 1 does not yield any meaningful result despite having the full-field response measurements. This case study will hence be limited to Method 2 only. Since only a single load source is present, in theory, a single sensor should be sufficient of estimating the load. However, it is needless to say that the number of sensors needed depends on the position of the sensors. If a sensor is positioned at a nodal point for a mode, the corresponding modal load will be poorly estimated.
Reduction in sensor information can be done through the selection matrix as follows:where will be the reduced measurement signal. The load scaling constant will consequently be
For example, if only the response is measured in DOF 1 and 5, the selection matrix becomes
This rewriting means that we can still identify loads in any of the six DOFs, although we only measure the response in a few selected DOFs. Figure 14 shows the estimated loads using only two sensors located in DOF 1 and DOF 5.
Figure 15 shows similar load estimate, but here, the response is measured at only the bottom node, in DOF 1. We see that even with one sensor, a reasonable estimate is obtained. Only the discontinuity at the impulse is off. Note that neither of the sensors in Figures 14 and 15 is located at a mode nodal point.
The precision of the load estimate is evaluated from different numbers of modes in the FRF matrix. The load is consequently estimated at a location where the response is not measured. This corresponds to a flawless modal expansion process [38, 39] or system identification from a set of roving sensors. In general, one may argue that having more modes in the FRF than sensors pose a challenge from a system identification point of view. However, we will not justify this here.
Data Availability
The data used in this study are generated from a simulated case study and can be made available upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the funding received from the Centre for Oil and Gas–DTU/Danish Hydrocarbon Research and Technology Centre (DHRTC).