Abstract
Knowledge on loads acting on a structure is important for analysis and design. There are many applications in which it is difficult to measure directly the dynamic loads acting on a component. In such situations, it may be possible to estimate the imposed loads through a measurement of the system output response. Load identification through output response measurement is an inverse problem that is not only ill-conditioned but, in general, leads to multiple solutions. Therefore, additional information such as the number and locations of the imposed loads must be provided ahead of time in order to allow for a unique solution. This paper focuses on cases where such information is not readily available and uses the concept of motion transmissibility for the identification of loads applied to a structure. The identification of loads through measurement of structural response at a finite number of optimally selected sensor locations is investigated. Optimum sensor locations are identified using the D-optimal design algorithm to provide the most precise load estimates based on acceleration measurements using accelerometers. Simulation results for multi-degree-of-freedom (MDOF) discrete and continuous systems are presented to illustrate the proposed technique. It is seen that the proposed approach is effective in determining not only the number of applied loads as well as their locations but also the magnitude of applied loads.
1. Introduction
The inverse-load identification problem has been investigated in several works. Most of the previous studies attempted to overcome ill-conditioning of the inverse problem by assuming that the load locations were known ahead of time and focused primarily on estimating magnitudes of applied loads. Some of the studies were based on a use of frequency or impulse response functions. In the frequency response function (FRF) based approaches, the magnitude of applied force is reconstructed using the FRF matrix and used with the system’s measured response (strain, acceleration, etc.) in the frequency domain. Due to ill-conditioning of the FRF matrix, any small errors in the measured response can cause large errors in the reconstructed load.
Recently, a review of different methods for overcoming ill-conditioning (in load identification) in the frequency domain was done by Hui et al. [1]. Some approaches used regularization (Tikhonov regularization, iterative regularization, and singular value decomposition), while others used dynamic programming in the time domain as proposed by Trujillo and Busby [2]. An experimental approach for load reconstruction has been presented by Doyle [3]. Genaro and Rade [4] used a sequential integration of measured acceleration signals to reconstruct the applied load magnitudes in the time-domain. Vishwakarma et al. [5] presented a Moore-Penrose pseudoinverse method which implements a least-square solution to reconstruct loads applied on a missile structure.
A study by Liu and Han [6] showed that prior knowledge of load locations can change the ill-posed problem into well-posed problem. In many situations, load locations are not always available which led researchers to work on a simultaneous determination of load magnitudes as well as prediction of load locations for different types of loads. Some researchers worked on impact load location and magnitude identification such as Khoo et al. [7]. It was shown experimentally that using the pseudoinverse method with a proper number and locations for measurement of system response leads to an accurate impact load identification. Other researchers worked on identification of harmonic loads. By measuring the transverse displacement, D’Cruz et al. [8] were able to determine not only the magnitude but also the location and the phase of a transverse harmonic point load applied on a viscoelastic plate. An application of Betti reciprocal theorem with a reference load case to estimate the magnitude and location for an assumed spatial load shape of the harmonic point load was presented by Moller [9]. Based on an optimization approach that can be implemented in time or frequency domain, Wang and Chiu [10] simultaneously predicted the amplitude and location of load applied on an arbitrary structure subjected to impact and harmonic loading.
Several studies done by Maia and his collaborators [11–16] have focused on an estimation of location, number, and magnitude of loads imposed on multi-degree-of-freedom systems using the concept of transmissibility. The reconstruction of loads is done in two successive phases. In the first phase, the location and number of applied loads are estimated by using a transmissibility model. In the second phase, the load vector is reconstructed by multiplying the inverse of the structural FRF matrix with the system’s measured response. This approach uses system response, such as accelerations, to predict the load magnitudes and locations. While this technique provides promising results, the question of sensor placement was not addressed and was left as user’s choice.
Practically, there are many locations on a structure where the accelerometers or strain gages can be mounted for measurement of the system response. Due to financial constraints and/or restrictions on potential sensor locations, the number and the locations of sensors are limited. In previous as well as recent works that use the concept of transmissibility for load prediction, the number of sensors used was addressed, but little attention was paid to their locations. The placements of sensors were left to the engineering experience or judgement of the user. According to Masroor and Zachary [17], the accuracy of load estimation is strongly influenced by location of sensors. They showed that a random placement of sensors increases problem ill-conditioning, whereas a proper selection of sensor locations decreases problem ill-conditioning and improves the accuracy of load estimation.
Through a numerical example, we show that an arbitrary selection of sensor locations (on a structure) is not the best strategy when it comes to prediction of load locations, especially for problems with multiple applied loads. To overcome this shortcoming, it is important that sensors be placed at optimum locations on the structure. In this work, a D-optimal algorithm is used to place the sensors optimally. Recently, Gupta and Dhingra [18, 19] used the D-optimal algorithm to identify accelerometer locations to estimate magnitudes of dynamic loads. However, since their approach assumes that the load locations are known in advance, their method is limited to certain applications. Another study by Augustine et al. [20] used the D-optimal algorithm to predict magnitude of moving loads wherein the locations of the moving load were assumed to be known ahead of time.
The limitation of the approaches presented by Gupta and Dhingra and Augustine et al. is that the number of loads and their locations are assumed to be known ahead of time. The only unknowns are load magnitudes. The novelty of this paper lies in the fact that using the concept of motion transmissibility, we can solve the load identification problem wherein all three load components: number of applied loads, load locations, and load magnitudes are unknown. We also provide an answer to the question of sensor placement for improved load prediction. This is especially important when multiple loads are applied. It is seen that the efficacy of load estimation is improved when sensors are placed at optimum locations. These optimum sensor locations are determined using the D-optimization technique.
Two examples dealing with numerical validation of the proposed approach are presented. The results of the examples indicate that the proposed approach is effective in determining number, locations, and magnitudes of loads applied to discrete as well as continuous structures.
2. Motion Transmissibility for MDOF Systems Based on FRF
Prior to 1998, the concept of transmissibility was largely limited to single-degree-of-freedom (SDOF) systems. In a multi-degree-of-freedom (MDOF) system, there are many input and output responses. Therefore, the transmissibility matrix is not unique for MDOF systems. Since 1998, the concept of transmissibility has been extended to MDOF by several researchers such as Riberio et al. [11] and Liu and Ewins [21]. Varoto and McConnell [22] discussed motion transmissibility concepts in the context of industrial applications and developed a matrix to characterize transmissibility of MDOF systems.
One of the approaches to determine the transmissibility of motion for MDOF systems is based on a use of FRF matrix which relates the dynamic displacement amplitudes D to the applied force amplitudes F. This matrix, also known as the receptance frequency response matrix, is given as follows:where is defined aswhere , , and are the stiffness, mass, and damping matrices, respectively, and can be generated from the finite element model of the structure. Rewriting equation (1) to recover the force vector in the frequency domain, we get
The transmissibility function is traditionally defined as the ratio of two different output spectra. For an MDOF system, it is better to divide the system coordinates into three groups as shown in Figure 1. Here, coordinates correspond to locations where the forces could be applied to the structure, whereas i coordinates are locations where the displacement responses Di are known or measured. The j coordinates are locations where the displacement responses Dj are unknown.
Rewriting equation (1) for i and j coordinates yields
From equation (5), using and substituting in equation (4) yieldswhere denotes the pseudoinverse of the FRF matrix, and the transmissibility matrix which relates both sets of displacements is defined as
The pseudoinverse of the submatrix can be obtained experimentally or analytically. The only required condition for the pseudoinverse to exist in equation (7) is the number of response measurements at i coordinates should be greater than or equal to the number of applied point loads at coordinates, i.e., .
To create the receptance matrix for the system under consideration, one can use the finite element method (FEM). The structures under consideration are modeled herein using ANSYS. The maximum number of modes returned from a finite element model is equal to model degrees of freedom, which can be a large number. Therefore, in many problems, numerical considerations make it impractical to retain all modes. Hence, a limited number of modes are retained that are “enough” to approximate the receptance matrix. In this work, the decision on the number modes retained depends on the cumulative mass fraction captured by retained modes. For reasonable accuracy, an adequate number of modes should be retained such that at least 90% of the cumulative mass fraction is captured by the retained modes.
3. Load Identification Using Motion Transmissibility with Arbitrary Placement of Sensors
The objective of the load identification problem is to estimate the locations and amplitudes of an unknown number of loads applied to a given structure using the system response data regardless of sensor placement on the structure. The solution process is divided into two phases as proposed by Lage et al. [12]. The first phase deals with the prediction of the number and the locations of the loads applied to the structure under consideration. This is done by searching for the right transmissibility matrix that matches the system dynamics (Neves and Maia [13]). The second phase involves generating the estimated load vector at the predicted load locations obtained in the first phase.
For the first phase, the number of the applied loads could be unknown along with their locations . Therefore, the search process for the transmissibility matrix transforms the dynamic responses into examining all possibilities until the predicted response matches the measured response . Based on the assumption made regarding the number of applied loads, various combinations of the test nodes are checked. For the case where it is assumed that only one load applied , the search process will start from the first node until the last node on the structure is traversed; the combinations of tested nodes will be . For the case with two applied loads , the combinations tested nodes will be . This approach can be extended to cover all possible combinations of load locations and number of applied loads .
The error between the measured dynamic response and the predicted dynamic response at each coordinate k is calculated as shown in equation (8). The error values are saved in an error vector as given in equation (9). They are analyzed later to find the combination with least associated error between the calculated response vector and the measured response vector :
The norm of is the assembled error for a specific combination of coordinates where forces F could possibly be located. To investigate all possible combinations for number of loads and their locations, the calculations are performed sequentially. The combination of load locations that gives the lowest error yields the number as well as the locations of the loads applied to the structure.
After determining the number and locations of the applied loads, the second phase involves solving the inverse problem by multiplying the pseudoinverse of the submatrix of FRF by the measured dynamic responses to reconstruct the load magnitudes:
An example is presented next that illustrates this idea. The results of this example will illustrate the need for optimum placement of sensors which is developed in Section 4.
3.1. Example 1: Discrete System with Two Applied Loads
To illustrate the load estimation using the concept of motion transmissibility discussed in Section 3, a 15-degree-of-freedom system shown in Figure 2 was analyzed with the following assumptions:(i)The system is undamped(ii)The first and the last mass are connected to fixed boundaries(iii)Masses are assigned arbitrary values starting from 20 kg for m1 to 160 kg for m15 in 10 kg increments(iv)Spring constants are assigned arbitrary values starting from 1 × 108 N/m for k1 to 8 × 108 N/m for k16 in increments of 0.5 × 108 N/m(v)Two sinusoidal forcing functions applied to masses m5 and m10: and
The first task is to determine the number and the locations of the applied loads using the localization method described earlier. Using modal analysis, it is found that the cumulative mass fraction (Irvine [23]) captured by the first five modes is 97% as shown in Table 1. Based on this observation, it was decided to retain five modes to reconstruct the receptance matrix . Sensor locations (i coordinates) are chosen to be uniformly distributed along system [3, 6, 9, 12, 15], while j coordinates [4, 13] are chosen randomly for applied load locations at coordinates [5, 10].
It is seen that the proposed approach can correctly find the number of applied loads, i.e., two applied loads. As can be seen from Figure 3, the accumulated errors have significantly low values for load combinations that correspond to two applied loads, which gives an accurate prediction for the number of loads applied. However, the locations of the two applied loads are not predicted accurately. The right combination of two applied loads at masses 5 and 10 is 70. Using five retained modes and a nonoptimal placement of sensors, the minimum error is seen at combination number 60 which corresponds to load location on masses 4 and 10.
Another attempt was made to improve the prediction of load locations by increasing the number of retained modes from 5 modes to all 15 modes. As can be seen in Figure 4, the minimum error occurs at combination number 72 which corresponds to the case when the loads are located at masses 6 and 10. It may be noted here that these minima have one load location correctly identified while the other one is near the actual location of node 5.
To examine the effect of sensor placement on prediction of load locations, several attempts were made for different sets of arbitrary locations for sensor positioning. The result for each set of sensor locations is shown in Table 2. Based on the results presented in Table 2, it can be seen that none of these attempts led to accurate load location prediction. An arbitrary selection of sensor locations is likely to get trapped at a local minimum of the error function and does not provide correct load locations. Since load identification is a two-phase sequential process, if the first phase does not yield accurate predictions for load locations, then the second phase is quite likely lead to inaccurate prediction for load magnitudes.
An important conclusion that can be drawn from this example is that to improve the localization approach using the transmissibility of motion, it is important to pick the locations of the sensors (i coordinates) carefully.
4. Load Identification Using Motion Transmissibility and Optimum Sensor Placement
As shown in the previous example, an arbitrary location of sensors to localize the applied loads may lead to inaccurate predictions of load locations. Typically, there are a large number of locations on the structure where the sensors (accelerometers) can be potentially mounted, and the precision with which the applied loads are estimated from measured sensor response is strongly influenced by the locations selected for sensor placement. A solution approach based on the construction of D-optimal designs is presented to determine the number and optimum locations of sensors that will provide the most precise load identification estimates.
4.1. D-Optimal Design for Sensor Locations
The objective of the load localization and load reconstruction problem using the concept of transmissibility of motion is that when a limited number of sensors are used (i coordinates), one should look for the best transmissibility that minimizes the error between predicted and actual response. An optimum placement of sensors can help accomplish this objective.
Using modal coordinates, the system displacements can be written aswhere is the modal matrix and denotes the mode participation factor (MPF). The modal matrix will be used in the D-optimal design to look for optimum sensor locations for measurement of output response, i.e., accelerations . Since , the least-squares approximation of is defined as
To understand the logic behind the D-optimal design in determining the optimum locations of the sensors, it is worthwhile to mention the basic idea behind this approach is to minimize the determinant of (ATA)−1 or equivalent that maximizes the determinant of the information matrix ATA of the design. For our problem, the system matrix A corresponds to the modal matrix . The solution approach involves maximizing determinant of . This is done using a MATLAB-based algorithm that sequentially adds and deletes points from a potential design by using a candidate set of points spaced over the region of interest. It is known that the acceleration vector is susceptible to measurement errors. Based on the statistical study of Masroor and Zachary [17], if the random errors in acceleration measurements are mutually independent and have the same standard deviation σ, then the variance-covariance matrix for the predicted load is given as
The matrix is known as the sensitivity of matrix . The precision of load estimates depends on the variance in the acceleration measurements and the conditioning of the sensitivity matrix. The accuracy of load estimates can be improved by improving the conditioning of matrix . Two factors that affect the sensitivity matrix are the number of sensors used and their locations on the structure. Therefore, choosing the optimum location and the suitable number of sensors can minimize the sensitivity of ; consequently, the variation in the load estimates will be minimized.
Kammer [24] and Atkinson et al. [25] have presented algorithms which help reduce the variance of matrix . The criterion of most relevance to the current application involves the maximization of the determinant of . Designs that maximize are called D-optimal designs. The D-optimal designs guarantee low variance among parameters to be estimated and low correlation between the estimated parameters.
A solution procedure is presented next that is used to provide the most precise estimates of the applied loads through optimum selection of the locations and the number of sensors (accelerometers) on the structure. The procedure, as shown in the flow chart in Figure 5, starts with determining the number of accelerometers to be used as . The general condition is that the number of accelerometers on i coordinates should be greater than or equal to the number of applied loads to be estimated , . Next, a candidate set is generated and examined to determine the i locations that provide the minimum variance in the load estimates. Based on the required number of accelerometers, the algorithm selects the optimum accelerometers from the candidate set which satisfy the condition stated above.
The goal is to generate an -point D-optimum design by selecting accelerometers’ locations from the candidate set such that the determinant of is maximized. From a randomly selected initial point design, the expansion algorithm adds a row of candidate point to the matrix such that the variance of predicted load is minimized. From the point design, optimal reduction of the expanded design is executed by deleting the row of candidate point that minimizes the variance of predicted load. This procedure of expanding and deleting candidate points continues until no further improvement in the objective function can be made. This is the central idea behind the exchange algorithm used herein.
The exchange algorithm implemented in MATLAB starts by choosing randomly different points from a given candidate set for a required number of loads that need to be estimated so the matrix of can be established. Next, a new candidate point is selected to be added as a row to expand the matrix and form so the determinant of is maximum. Out of the rows in matrix , a row is deleted to construct such that the determinant of is maximum. This process of adding and deleting rows continues until the determinant value of cannot be improved further. The final result of the D-optimization, , is the matrix which provides the information on the optimum accelerometer locations that will provide minimum variance in the load estimates.
As explained above, the process of adding and deleting rows from matrix requires calculating the determinant at each step. This calculation is numerically expensive and should be minimized. An alternate method for calculating the determinant = from that of when a row is added to the matrix iswhere [+] denotes addition and is substituted by subtraction [−] when row is deleted from . To be able to use equation (14), needs to be updated also. This updating when row is added to the matrix is given bywhere [−] denotes subtraction and is substituted by addition when deleting a row from .
With the optimum locations of the accelerometers determined , the accelerometers are mounted at those locations before the application of the unknown loads. The measured accelerations vector , along with the optimum , is then used to estimate the unknown forces as shown equation (10).
Two examples are presented next that illustrate numerical validation of the approach presented above. The examples also illustrate the effectiveness of using optimum sensor locations to identify the loads applied to a structure.
4.2. Discrete System Example Revisited: Load Identification Using Optimum Sensor Locations
The numerical example dealing with 15-DOF spring-mass system described earlier is revisited to identify imposed loads using D-optimal design for optimum placement of sensors. Using the exchange algorithm outlined in Figure 5, optimum locations for the sensors were determined as i = [4, 6, 8, 11, 15] while j = [3, 9] are chosen randomly for same load locations [5, 10]. The result in Figure 6 shows that by using the optimum accelerometer locations, the minimum accumulated error occurs at the right combination number (70) for two applied loads at masses 5 and 10. An accurate prediction of load locations will lead to accurate predictions of the load magnitudes as shown in Figures 7 and 8.
To simulate an experimental situation where the accelerations are measured experimentally and measurement errors may be present, each element of response measurement in was corrupted with normally distributed random errors with zero mean and standard deviation of 10% of its value.
Using the algorithm described above, five optimum locations for accelerometers at i coordinates are found to be [2, 6, 8, 12, 14], and j coordinates [3, 9] were chosen arbitrarily for same load locations [5, 10]. The results in Figure 9 show that by using the optimal locations for accelerometers and with errors present in response measurements, the minimum accumulated error occurs at the right combination number (70) for two applied loads at masses 5 and 10. It may be noted that due to the presence of measurement errors in accelerometer readings, the absolute values of the accumulated errors have increased. As shown in Figures 10 and 11, it can be seen that the applied loads are recovered accurately despite the presence of measurement errors.
Based on a close agreement between applied and predicted loads, the results of this example indicate that the proposed approach was in this case effective in not only determining load magnitudes but also unknown load locations. An example dealing with an application of proposed approach to a continuous system is presented next.
4.3. Example 2: Load Identification on a 3D Cantilever Beam
The numerical example discussed previously dealt with a discrete system. Next, a continuous system is considered where dynamic load prediction is applied to a 0.25 m long, 0.05 m wide, and 0.005 m thick cantilevered steel beam (see Figure 12). Without loss of generality, the system is assumed to be undamped. A finite element model of the beam was developed using ANSYS and meshed with SHELL181 elements. The beam has 36 nodes and each node has six degrees of freedom. Six of these nodes are completely constrained; so, the structure has 30 unconstrained nodes with 180 degrees of freedom.
At the free end of the beam, a vertical load is applied on node 19 described as with another vertical load applied on node 24 given as .
As shown earlier for the discrete system, the solution approach will involve a prediction of the locations of the applied loads followed by a reconstruction of the load magnitudes by using the motion transmissibility and optimum locations for the sensors. A modal analysis was performed on the FE model of the beam to obtain the modal matrix for the structure. For this 180-DOF example, the modal matrix is . The mass and stiffness matrices were obtained using finite element method. ANSYS provides data for and matrices in the Harwell-Boeing file format. A program was written in MATALB to convert them into the matrix format suitable for current application.
As discussed earlier, a limited subset of modes is retained to reconstruct the applied loads. The retained modes should capture at least 90% of the cumulative mass fraction.
If only m modes are retained to reconstruct system response, the condensed modal matrix is an matrix. If the direction of the applied loads is known a priori, then as a first step, it may be adequate to construct the reduced modal matrix such that it has only the modes in the same direction as the applied loads. For this example, the reduced modal matrix will have thirty normal modes in Y direction, . Therefore, the candidate modal matrix will be , and this will be the input for D-optimal program.
Following the D-optimal design algorithm described previously, the candidate set is searched to determine its optimum subset . After is found, the optimum accelerometer locations are determined. The accelerometers are mounted at the identified optimum locations and the acceleration is measured, which can be then be successively integrated numerically to obtain and . By using the finite element model in ANSYS, the displacement vector can be found directly from ANSYS at the optimal locations. From , one can use MATLAB program to get the responses in the frequency domain which represent the displacement vector at optimal i coordinate and will be used later in the force reconstruction step as shown in equation (10). The optimal (i) coordinates for the cantilevered beam are found to be [3, 4, 10, 12, 16, 33].
To find the locations of the applied forces, both measured and predicted displacement vectors at j coordinates should be known to look for the minimum error as in equation (8). For the measured responses at j coordinates, one can arbitrarily pick any locations. For the example considered here, j coordinates are assumed to be same i coordinates.
For the predicted response at j coordinates, all possibilities were explored until the calculated response matched the measured ones. This method was implemented in MATLAB. The algorithm scrolls through possible combinations of applied force locations. For each combination, it calculates the associated error between the calculated vector and the measured responses vector; this is done over the range of frequencies defined by the user. For each combination, the calculated error is saved in an error vector and plotted as shown in Figure 13. For this example, the total number of combinations explored for the case of two applied loads is 465. The applied load locations at nodes 19 and 24 correspond to combination number 364. It can be seen from Figure 13 that there is a minimum value at this combination number which corresponds to the load location being predicted correctly.
It is worth mentioning that in Figure 13, there are three other minima which belong to the following combination numbers: (i) 276 which corresponds to the case when two loads are applied on nodes 14 and 15, (ii) combination number 284 which corresponds to loads on nodes 14 and 24, and (iii) combination number 298 and which corresponds to loads being applied on nodes 15 and 19. It can be seen from Figure 12 that node 19 lies above node 14 while node 24 lies above node 15. Since both nodal pairs share the same applied load location, the algorithm is likely to pick one node from each of the two pairs.
During the next step, the load magnitudes were reconstructed using equation (10) and transformed into time domain using inverse Fourier transform. The reconstructed loads are plotted along with applied loads as shown in Figures 14 and 15. It can be seen from both figures that the load trends are recovered with reasonable accuracy.
5. Conclusions and Future Work
In this paper, load identification (load location as well as magnitude) by using the concept of motion transmissibility has been examined for two different multi-degree-of-freedom systems: a discrete system and a continuous system. Based on the results presented, it is shown that to improve the accuracy of the load location prediction problem, the placement of sensors at correct locations is important. Using optimum locations of accelerometers as determined by the D-optimal algorithm is able to improve the identification for the unknown loads especially when multiple loads are applied and when the error function has multiple local minima. In the present study, it has been verified numerically that even in the presence of simulated measurement errors, the proposed method was able to achieve promising results. Future work on this approach will include an experimental validation of this method. It was seen that increasing the number of retained modes to reconstruct the response improved the accuracy of the obtained load identification results. Practically, there are limitations on the number of modes whose MPF can be estimated from accelerometer measurements. Further research on modal condensation shall be conducted to overcome this limitation and so that the accuracy of load identification results can be improved further.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.