Abstract
This paper develops a technique for identifying dynamic loads acting on a structure based on impulse response of the structure, also referred to as the system Markov parameters, and structure response measured at optimally placed sensors on the structure. Inverse Markov parameters are computed from the forward Markov parameters using a linear prediction algorithm and have the roles of input and output reversed. The applied loads are then reconstructed by convolving the inverse Markov parameters with the system response to the loads measured at optimal locations on the structure. The structure essentially acts as its own load transducer. It has been noted that the computation of inverse Markov parameters, like most other inverse problems, is ill-conditioned which causes their convolution with the measured response to become quite sensitive to errors in system modeling and response measurements. The computation of inverse Markov parameters (and thereby the quality of load estimates) depends on the locations of sensors on the structure. To ensure that the computation of inverse Markov parameters is well-conditioned, a solution approach, based on the construction of -optimal designs, is presented to determine the optimal sensor locations such that precise load estimates are obtained.
1. Introduction
It is often desired to have knowledge of loads acting on a structure early in design process so that effective analysis and design optimization, ensuring the structural integrity of the structure, may be performed. Precise identification of the loads is imperative to having high level of confidence in numerical simulation, such as finite element analysis, which precludes the need for expensive and time consuming experimental testing. In some cases, the loads acting on a structure may be measured by introducing load transducers (load cells) between the structure and the applied loads. This method of load measurement has some disadvantages. Introduction of load transducers may change the system dynamic characteristics and the system may not behave in the same fashion as it would if the load transducers were not present. This leads to inaccuracies in the load estimation. In some applications, the input locations may not be accessible to insert a load transducer for the measurement of loads being transmitted to the structure.
The limitations posed by external load transducers may be overcome by transforming the structure into a self-transducer, where instead of measuring the loads directly, response of the structure to the unknown applied loads is measured. Physically measurable quantities such as displacements, accelerations, and strains that depend on the loads can be used as response. Some kind of relationship between the measured quantity and the loads to be estimated is then utilized to estimate the applied loads. The instrumented structure effectively becomes a load transducer. This class of problems is known as the “inverse problem.”
There are several challenges associated with solving the inverse problem. The major challenge is that the inverse problems tend to be highly ill-conditioned. The difficulties posed by the inverse problems pertaining to load identification are documented by Stevens [1]. Ewins [2] used displacement transducers to identify input loads and found that very small variations (noise) in the response measurement cause large errors in the force estimation. Load identification problem in frequency domain was studied by Desanghere and Snoeys [3] where they attributed the reason for ill-conditioning to a few dominant elements in the frequency response function (FRF) matrix. Hillary and Ewins [4] and Okubo et al. [5] investigated the influence of noise in the measured response and noted that the inversion process is highly ill-conditioned.
Another technique was proposed by Busby and Trujillo [6] where the load estimation problem is treated as minimization problem of error between the measured structural response and the response predicted from model, which is then solved using dynamic programming. This technique was extended by Hollandsworth and Busby [7] by applying it to actual experimental measurements. Application of the technique is limited since the amount of computation increases dramatically with the model order.
Genaro and Rade [8] suggested a load reconstruction method from acceleration response and its successive integration. Carne et al. [9] developed a technique for load identification known as the Sum of Weighted Accelerations Technique (SWAT). An optimization problem for load recovery was formulated by Hashemi and Kargarnovin [10] where the objective function is the difference between analytical and measured responses and the decision variables are the location and magnitude of the applied force. Lu and Law [11] proposed a load estimation technique based on sensitivity of structural responses. A gradient-based model updating method based on dynamic response sensitivity is utilized to identify the force and physical parameters. Boukria et al. [12] applied the FRF technique to identify the location and magnitude of impact force on a circular plate.
Lee [13] presented an impact load recovery technique on thick plates utilizing Green’s function and acoustic waveform due to impact load, which are theoretically obtained based on plate theory. Load identification method on a nonconservative nonlinear system which has a nonlinear damping effect was studied by Jang [14]. Ahmari and Yang [15] proposed an inverse analysis technique based on mathematical forward solving model for a simply supported plate and Particle Swarm Optimization (PSO) to identify impact location and impact load time history. Khoo et al. [16] developed an impact force identification technique utilizing Operating Deflection Shape (ODS) analysis, FRF measurement, and pseudoinverse method. They concluded that the accuracy of identified loads depends on sensor locations and improves with increased number of measured response locations. Liu et al. [17] proposed load reconstruction technique in time domain based on shape function method of moving least-square fitting and regularization.
The reason for ill-conditioning of the load prediction problem in frequency domain was investigated by Starkey and Merrill [18] where they attributed it to near-singularity of the FRF matrix. Hansen and Starkey [19] furthered the investigation on the nature of ill-conditioning in modal domain. They studied the effect of locations of accelerometers on a steel beam on the condition number of the modal matrix and concluded that the condition number of the modal matrix can be improved through proper placement of accelerometers on the beam and proper selection of modes included in the analysis.
Kammer [20] presented a method that utilizes a set of inverse system Markov parameters estimated from forward system Markov parameters using a linear predictive scheme. This computation is ill-conditioned; therefore, a regularization technique is employed to stabilize the computation. Steltzner and Kammer [21] suggested a technique for input force estimation using an Inverse Structural Filter (ISF) that processes the structural response data and returns an estimate of the input forces. In the aforementioned works, the authors made no mention of any strategy to optimally place sensors on the structure that would result in minimizing the effect of ill-conditioning associated with the inversion. A number of these works [8, 9, 20] assume that the accelerometers are collocated with the forces, which may not always be feasible.
A study on the effect of strain gage locations and angular orientations on the quality of load estimates was conducted by Masroor and Zachary [22]. Through formulation of a sensitivity parameter, they argued that proper selection of gage locations improves the condition of the inverse problem. Since analysis based on all possible gage location combinations is not feasible, they selected a few groups of gages based on judgment for the analysis. This approach may still result in poor load estimates because the sensors may be placed at a location where they have relatively low sensitivity to the applied loads.
In view of the above discussion, it can be concluded that although a lot of researches have been dedicated towards identifying novel techniques for load estimation, the idea of determining optimal sensor locations, so as to minimize the effect of modeling and measurement errors, has received little attention. It is recognized [16, 19, 22] that sensor location has significant influence on the overall quality of results. To overcome the abovementioned shortcomings, this work expands upon that by Kammer [20] and Steltzner and Kammer [21] by proposing an approach for estimating time-varying loads acting on a structure by measuring structure response at finite number of optimum locations on the structure. The accuracy of the load estimates is dependent on the locations of the sensors. The novel approach, based on -optimal design, results in increased accuracy in the dynamic load estimation. Two examples dealing with numerical validation of the proposed approach are presented to illustrate the effectiveness of the proposed technique.
2. Theoretical Development
2.1. Markov Parameter Representation
Any continuous system may be spatially discretized as an n degrees of freedom (DOF) system, whereby its equilibrium equation of motion can be written in matrix form as For simple systems, , , and may be obtained by writing the system equations of motion; for complex systems, they can be generated from the finite element model of the structure.
Another way to characterize the input-output behavior of a structural system is by writing the second-order (1) as a first-order equation or in continuous time-invariant state-space form as where is the system matrix and is the input matrix for continuous case given by
All physical phenomena fundamentally exist in continuous time. The experimentally measured response data, however, is available only at discrete time instants. This calls for the need to transform the continuous time-invariant state-space model into a discrete time-invariant state-space model. The transformation takes the following form: where is the subscript over the discretized time and and are related to their continuous case counterparts as Equation (4) is the discrete linear time-invariant state-space representation of the structural system. The corresponding system output is given by Given that unit impulse load is applied to the system at time , that is, , and assuming zero initial conditions, the impulse response at various time points is given by Given the impulse response, the response at any time is given by the convolution of the input force with the impulse response as where the matrices are called forward Markov parameters [20] and summarized as
The forward Markov parameters represent the response of the discrete system to applied unit impulse and thus contain the dynamic properties of the system. They can be obtained either analytically from the discretized model of the structural system or experimentally by measuring the output of the system due to a known input, computing the corresponding frequency response function, and then taking its inverse discrete Fourier transform. Equation (8) is the Markov parameter representation of a structural system.
2.2. Problem Formulation
In the inverse problem of interest here, the objective is to estimate the input forces , given the forward Markov parameters and system responses.
From (8), the system response at is given by Assuming the case where the number of sensors used is at least equal to the number of loads, the least-squares estimate of the loads at is given as Similarly, at , the system response is given by which can again be solved for the loads as By induction, the loads at time ti can be verified to be
Next, define a discrete inverse system similar to (4) and (6), where the roles of input forces and system response are reversed, as where the inverse matrices , , , and are related to the corresponding forward matrices by Similar to (8), given the system response, the input forces at any time are given by the convolution relation as where the matrices are known as the inverse Markov parameters [20]. The inverse Markov parameters, similar to the forward Markov parameters, contain the dynamic properties of the inverse system. On comparison of expansions of (14) and (17), it may be ascertained that the inverse Markov parameters are related to the forward Markov parameters by a linear predictive equation given by
There exist cases where is zero matrix. In certain other nonminimum phase structural systems, is rank deficient and the least-squares inverse in (18) does not exist. This renders the causal summation in (17) undefined. To deal with this limitation, the system output (6) is stepped forward in time as [21] The inverse system associated with (19) is noncausal; that is, the input force estimates at current time become a function of the structural response at future times. In general, if the first z Markov parameters in (7) are zero, that is, , the noncausal z-lead inverse system can be constructed as and the corresponding force estimates are given by It must be noted in (21) that the input force estimates at current time ti are dependent on the structural response at future times .
The determination of the inverse Markov parameters from the forward Markov parameters using (20) is computationally expensive; however, it must be noted that the computation of the inverse Markov parameters needs to be performed only once for any particular system. Once the inverse Markov parameters are computed, the input forces can be predicted from the response using (21). This technique, however, suffers from a similar limitation as most of the other inverse problems—ill-conditioning. The inaccuracies in system modeling translate to errors in forward Markov parameters. The computation of inverse Markov parameters from erroneous forward Markov parameters is ill-conditioned, and the computed inverse Markov parameters become unbounded. Thus, the convolution of the inverse Markov parameters with the noisy structural response given by (21) is not always a converging sum. The sum becomes numerically unstable due to the intrinsic ill-conditioning of the inverse problem identified by (20).
Conditioning of the erroneous forward Markov parameters depends on the locations of sensors on the structure; therefore, the input loads can be estimated precisely from measured structural response by strategic placement of sensors on the structure. Since it is not feasible to place sensors at all the DOFs of a structure, a methodology is presented next to overcome the identified shortcoming.
2.3. Candidate Set
Optimal selection of sensor locations on the structure requires a candidate set to be defined. It can be inferred from the load identification step in (21) that the quality of load estimates depends directly on how accurately the inverse Markov parameters are computed. It can further be assessed from (20) that the th inverse Markov parameter is computed by the least-squares inversion of the th Markov parameter . All the other higher order inverse Markov parameters depend on . Therefore, the accuracy in the computation of from dictates the accuracy of all the other higher order inverse Markov parameters. Each row in corresponds to a unique DOF of the system output. Typically, there are a large number of locations on a structure where sensors can potentially be mounted. These locations may exclude inaccessible locations such as the regions of load application. To form a candidate set for optimization, the DOFs (rows) corresponding to the inaccessible locations for sensor placement are ignored from ; the subset so obtained is called the candidate set .
To obtain a meaningful solution to (20), an optimal subset of the candidate set needs to be identified, which in turn will lead to an increased accuracy in the estimation of . As more sensors are used, the additional information helps to obtain a more precise estimate of , but practical and financial constraints place limitations on the number of sensors that can be used. If the number of forces to be estimated is , then in order to minimize the error in estimates, the number of sensors a must satisfy the criterion . Each sensor location, referred to as a candidate point, provides a potential row for inclusion in . needs to be such a subset of that provides the most precise estimates of . An approach to extract from is described next.
2.4. -Optimal Design
If the errors in response measurements are statistically independent and the standard deviation of each of them isσ, then the covariance matrix for estimates is given by [22] The matrix is known as the sensitivity of . For a given variance in response measurements, minimization of the sensitivity of leads to an increased precision in estimates. The sensitivity of is a function of the number and locations of the sensors on the structure. Therefore, optimum selection of the number and locations of the sensors on the structure so as to minimize the sensitivity of leads to the minimization of variation in the estimates.
For a given number of sensors a, the candidate set is searched to determine a sensor locations that provide the least variance in the load estimates. One of the criteria followed involves the maximization of the determinant of . Designs that maximize are called -optimal designs [23–25], where D denotes determinant. In order to construct an a-point -optimal design, the a sensor locations that maximize must be selected from the candidate set . To select the a-point -optimal design, sequential exchange algorithms [24] based on the principles of optimal augmentation and reduction of an existing design are used herein.
The basic idea behind the sequential exchange algorithm is as follows. Given the candidate set and the number of sensors a, the first step is to randomly select a distinct candidate points from the candidate set to initialize . Out of the remaining candidate set, a candidate point is then selected and the corresponding row is augmented to to form such that is maximum. Next, out of the rows in , a row is deleted to arrive at such that is maximum. This process of adding and deleting rows continues until there is no further improvement in the value of . The final so obtained is the -optimal design whose row indices provide the information on the optimum sensor locations.
As the candidate set gets bigger and bigger, increasingly large number of determinants and matrix inverses need to be calculated. A very naïve way to compute the determinant is to obtain and then calculate the determinant . This approach of determinant calculation is computationally very expensive. An alternate formula for computing the determinant from that of when a row is augmented to is [25] where denotes addition. The is replaced by subtraction for the case when a row is deleted from . In order to be able to use (23), can be maintained and updated as the row is augmented to by where denotes subtraction and is replaced by addition for the case when a row is deleted from .
Once and, in turn, optimum sensor locations are determined, sensors are mounted at the identified optimum locations and response is measured. Equations (20) and (21) are then utilized to compute the input forces . It is to be noted that in this approach it is not necessary that response measurements be taken at locations where forces are applied.
2.5. Error Quantification
The percent root mean square (rms) error in the recovered loads with respect to the actual applied loads can be calculated as The above equation can be used to quantify errors in the load estimates.
3. Example 1: 15-DOF Spring-Mass-Damper System
The dynamic load estimation method from acceleration measurements is illustrated with the help of a numerical simulation consisting of a 15 DOF spring-mass-damper system shown in Figure 1. Masses and are connected to fixed boundary. A periodic forcing function is applied to mass . The task is to determine the optimum accelerometer locations and reconstruct the input force based on the acceleration time response at those locations.
A total of 2 accelerometers were used to measure accelerations of 2 masses. The , , and matrices were obtained by writing the equations of motion for the system. Discrete time-invariant state-space form of the system was obtained using (3)–(6). Forward Markov parameters were computed from the state-space model using (9). All numerical computations were performed in a MATLAB programming environment.
All DOFs, except the DOF where the load was applied, were selected to be the locations where accelerometers can potentially be mounted; that is, the DOF corresponding to the applied load did not form a part of the candidate set . was subjected to the -optimal design algorithm described in Section 2.4 to obtain whereby the optimum accelerometer locations were found to be at masses and . Having obtained the optimum accelerometer locations, the inverse Markov parameters were computed using (20). In absence of any experimental data, the acceleration time response at the optimum locations was obtained by solving the ordinary differential equations numerically. Using the acceleration data computed at the optimum accelerometer locations, the input force was recovered using (21). The applied and the recovered loads are plotted in Figure 2. The rms error using (25) was calculated to be 0.7%. It can be inferred from the plot that the applied load is recovered accurately.
4. Example 2: Overhead Crane Girder
A more general numerical example is presented here where a vertical dynamic load acting on an overhead crane girder, through trolley wheels, needs to be estimated. For the sake of illustration and simplicity, the number of loads is assumed to be one and assumed to be applied at the girder mid-span. The problem is to determine the optimum accelerometer locations and reconstruct the input force based on the acceleration time response at those locations.
A finite element model of the girder was developed in ANSYS using BEAM188 elements. The finite element model of the girder along with the applied load and boundary conditions is shown in Figure 3. The model consisted of 50 beam elements and 49 unconstrained nodes with 6 degrees of freedom per node; that is, the total number of degrees of freedom of the model was 294.
The and matrices were obtained using finite element method. ANSYS provides data for and matrices in the Harwell-Boeing file format. A routine was written in MATLAB to convert them into the matrix format suitable for current application. All further calculations were performed in MATLAB. Discrete time-invariant state-space form of the system was obtained using (3)–(6). Forward Markov parameters were computed from the state-space model using (9). The number of accelerometers to be used was arbitrarily assumed to be equal to 3. All DOFs, except the DOF where the load was applied, were selected to be the locations where accelerometers can potentially be mounted; that is, the DOF corresponding to the applied load did not form a part of the candidate set . was subjected to the -optimal design algorithm described in Section 2.4 to obtain . The optimum accelerometer locations are also depicted in Figure 3. It is observed that the algorithm predicts the optimum sensor locations to be as close to the loads as possible. This should not be considered as a limitation to the application of the proposed technique since if certain locations around the force application points are not available for sensor placement, they can initially be excluded from the candidate set . Furthermore, if the sensor positions seem to be too congested mutually, additional criteria may be instructed such as if a particular spot is chosen as a potential sensor location, then certain area around that location may be excluded from the candidate set.
Having obtained the optimum accelerometer locations, the inverse Markov parameters were computed using (20). In absence of any experimental data, the acceleration time response at the optimum locations was obtained by solving the ordinary differential equations numerically. Using the acceleration data computed at the optimum accelerometer locations, the input force was recovered using (21). The applied and the recovered loads are plotted in Figure 4. The rms error using (25) was calculated to be 0.1%.
In order to simulate a more realistic scenario where accelerations are measured experimentally, the acceleration data was corrupted with normally distributed random errors with zero mean and standard deviation of 10% of its value. The applied and recovered loads, with errors in acceleration measurements, are plotted in Figure 5. The rms error using (25) was calculated to be 2.3%.
Next, a different nature of force was considered where a half-sine wave loading was applied. Both exact acceleration data and acceleration data corrupted with normally distributed random errors with zero mean and standard deviation of 10% of its value were considered. The applied and recovered loads, without and with errors in acceleration measurements, are plotted in Figures 6 and 7, respectively. The rms error using (25) was calculated to be 1.3% for the measurement involving acceleration error. It can be inferred from the plots that the proposed approach is able to successfully recover the applied load precisely even when realistic measurement errors are present in the structural response.
Based on the results for both examples, it can be seen that, for discrete as well as continuous systems, the proposed approach is able to accurately estimate the loads acting on a component by measuring the acceleration response at a finite number of optimum locations.
5. Conclusions
A computational technique is presented that allows for indirect measurement of dynamic loads acting on a structure such that the structure behaves as its own load transducer. This is achieved by placing finite number of sensors on the structure. The loads exciting a structure are estimated by convolving the structural response with the inverse Markov parameters. The inverse Markov parameters, in turn, are computed from the forward Markov parameters using a linear prediction algorithm. The forward Markov parameters represent the response of the system to applied unit impulse and thus contain the dynamic properties of the system. They can be obtained analytically as well as experimentally.
The computation of the inverse Markov parameters from forward Markov parameters, like all inverse problems, suffers from ill-conditioning. The computed inverse Markov parameters are not always bounded, thus rendering their convolution with the structural response to diverge. The accuracy of the inverse Markov parameters (and thereby the quality of input load estimates) depends on the locations of sensors on the structure. To improve the precision of load estimates, the technique of -optimal design is utilized to determine the optimum sensor locations. It is observed that the algorithm predicts the optimum sensor locations to be as close to the loads as possible. This should not be considered as a limitation to the application of the proposed technique since if certain locations around the force application points are not available for sensor placement, they can initially be excluded from the candidate set. Furthermore, if the sensor positions seem to be too congested mutually, additional criteria may be instructed such as if a particular spot is chosen as a potential sensor location, then certain area around that location may be excluded from the candidate set. The loads recovered based on accelerations measured from optimally placed accelerometers on the structure are observed to be in excellent agreement with the applied loads.
Nomenclature
| : | Zero matrix |
| : | Derivative with respect to time |
| : | System matrix for continuous case |
| : | System matrix for discrete case |
| : | Inverse system matrix for discrete case |
| : | Number of sensors |
| : | Input matrix for continuous case |
| : | Input matrix for discrete case |
| : | Inverse input matrix for discrete case |
| : | Damping matrix |
| : | Output matrix for discrete case |
| : | Inverse output matrix for discrete case |
| : | Feedforward matrix for discrete case |
| : | Inverse feedforward matrix for discrete case |
| : | Load vector |
| : | Applied load |
| : | Reconstructed load |
| : | th forward Markov parameters matrix |
| : | Candidate set |
| : | Optimum subset of determined by -optimal design |
| : | th inverse Markov parameters matrix |
| : | Identity matrix |
| : | Stiffness matrix |
| : | Mass matrix |
| : | Number of degrees of freedom of the system |
| : | Number of applied forces |
| : | Time |
| : | Discretized time |
| : | Time step |
| , : | State variables vector |
| : | Displacement vector |
| : | System output |
| : | Unit impulse load function |
| : | Percent root mean square (rms) error. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.