# Load Reconstruction Technique UsingD-Optimal Design and Markov Parameters

Authors: Deepak K. Gupta; Anoop K. Dhingra
Journal: Shock and Vibration
 (2015)

---

## 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 ofD-optimal designs, is presented to determine the optimal sensor locations such that precise load estimates are obtained.


---

## Body

## 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 D-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 ann degrees of freedom (DOF) system, whereby its equilibrium equation of motion can be written in matrix form as
(1)Mx¨t+Cx˙t+Kxt=ft.
For simple systems, M, C, and K 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
(2)u˙t=Acut+Bcft,
where Ac is the system matrix and Bc is the input matrix for continuous case given by
(3)Ac=0I-M-1K-M-1C,Bc=0M-1.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:(4)uti+1=Aduti+Bdfti,
where ti is the subscript over the discretized time and Ad and Bd are related to their continuous case counterparts as
(5)Ad=e[Ac]Δt,Bd=Ac-1Ad-IBc.
Equation (4) is the discrete linear time-invariant state-space representation of the structural system. The corresponding system output is given by
(6)yti=Cduti+Ddfti.
Given that unit impulse load is applied to the system at time t=0, that is, δ0=1, and assuming zero initial conditions, the impulse response at various time points is given by
(7)H0=Cdu0+Ddδ0=Ddu1=Adu0+Bdδ0=BdH1=Cdu1+Ddδ1=CdBdu2=Adu1+Bdδ1=AdBdH2=Cdu2+Ddδ2=CdAdBdu3=Adu2+Bdδ2=Ad2BdH3=Cdu3+Ddδ3=CdAd2Bd⋮Hi=CdAdi-1Bd.
Given the impulse response, the response at any time is given by the convolution of the input force with the impulse response as
(8)yti=∑i=0tiHifti-i,
where the matrices Hi are called forward Markov parameters [20] and summarized as
(9)H0=DdHi=CdAdi-1Bd,i=1,2,3….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 forcesf, given the forward Markov parameters and system responses.From (8), the system response at ti=0 is given by
(10)y0=H0f0.
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 ti=0 is given as
(11)f0=H0TH0-1H0Ty0.
Similarly, at ti=1, the system response is given by
(12)y1=H0f1+H1f0
which can again be solved for the loads as
(13)f1=H0TH0-1H0Ty1-H1f0.
By induction, the loads at timeti can be verified to be
(14)fti=H0TH0-1H0Tyti-∑i=1tiHtifti-i.Next, define a discrete inverse system similar to (4) and (6), where the roles of input forces f and system response y are reversed, as
(15)uti+1=AdIuti+BdIyti,fti=CdIuti+DdIyti,
where the inverse matrices AdI, BdI, CdI, and DdI are related to the corresponding forward matrices by
(16)AdI=Ad-BdDdTDd-1DdTCd,BdI=BdDdTDd-1DdT,CdI=-DdTDd-1DdTCd,DdI=DdTDd-1DdT.
Similar to (8), given the system response, the input forces at any time are given by the convolution relation as
(17)fti=∑i=0tihiyti-i,
where the matrices hi 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 hi are related to the forward Markov parameters Hi by a linear predictive equation given by
(18)h0=H0TH0-1H0T,hk=-h0∑i=1kHihk-i.There exist cases whereH0 is zero matrix. In certain other nonminimum phase structural systems, H0 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]
(19)yti=Cduti+CdBdfti.
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 firstz Markov parameters in (7) are zero, that is, Dd=CdBd=CdAdBd=⋯=CdAdz-2Bd=0, the noncausalz-lead inverse system can be constructed as
(20)hz=HzTHz-1HzT,hk+z=-hz∑i=1kHi+zhk-i+z
and the corresponding force estimates are given by
(21)fti=∑i=0tihi+zyti-i+z.
It must be noted in (21) that the input force estimates f at current timeti are dependent on the structural response y at future times ti-i+z.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 h with the noisy structural response y 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 loadsf can be estimated precisely from measured structural response y 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 zth inverse Markov parameter hz is computed by the least-squares inversion of the zth Markov parameter Hz. All the other higher order inverse Markov parameters depend on hz. Therefore, the accuracy in the computation of hz from Hz dictates the accuracy of all the other higher order inverse Markov parameters. Each row in Hz 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 Hz; the subset so obtained is called the candidate set Hz_cs.To obtain a meaningful solution to (20), an optimal subset Hz_opt of the candidate set Hz_cs needs to be identified, which in turn will lead to an increased accuracy in the estimation of f. As more sensors are used, the additional information helps to obtain a more precise estimate of f, 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 nf, then in order to minimize the error in f estimates, the number of sensorsa must satisfy the criterion a≥nf. Each sensor location, referred to as a candidate point, provides a potential row for inclusion in Hz_opt. Hz_opt needs to be such a subset of Hz_cs that provides the most precise estimates of fti. An approach to extract Hz_opt from Hz_cs is described next.


### 2.4.D-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 fti estimates is given by [22]
(22)varfti=σ2HzTHz-1.
The matrix HzTHz-1 is known as the sensitivity of Hz. For a given variance σ2 in response measurements, minimization of the sensitivity of Hz leads to an increased precision in fti estimates. The sensitivity of Hz 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 Hz leads to the minimization of variation in the fti estimates.For a given number of sensorsa, the candidate set Hz_cs is searched to determinea sensor locations that provide the least variance in the load estimates. One of the criteria followed involves the maximization of the determinant of HzTHz,HzTHz. Designs that maximize HzTHz are called D-optimal designs [23–25], whereD denotes determinant. In order to construct ana-point D-optimal design, thea sensor locations that maximize HzTHz must be selected from the candidate set Hz_cs. To select thea-point D-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 setHz_cs and the number of sensorsa, the first step is to randomly selecta distinct candidate points from the candidate set to initialize Hz. Out of the remaining candidate set, a candidate point is then selected and the corresponding row is augmented to Hz to form Hz+ such that HzT+Hz+ is maximum. Next, out of the a+1 rows in Hz+, a row is deleted to arrive at Hz- such that HzT-Hz- is maximum. This process of adding and deleting rows continues until there is no further improvement in the value of HzTHz. The final Hz so obtained is the D-optimal design Hz_opt 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 obtainM=HzTHz and then calculate the determinant M. This approach of determinant calculation is computationally very expensive. An alternate formula for computing the determinant M+=HzT+Hz+ from that of M when a row yT is augmented to Hz is [25]
(23)M+=M1+yTM-1y,
where + denotes addition. The + is replaced by subtraction for the case when a row yT is deleted from Hz+. In order to be able to use (23), M-1 can be maintained and updated as the row yT is augmented to Hz by
(24)M+-1=M-1-M-1yM-1yT1+yTM-1y,
where - denotes subtraction and is replaced by addition for the case when a row yT is deleted from Hz+.OnceHz_opt 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 f. 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(25)ε%rms=∑fapp-frec2∑fapp2×100%.
The above equation can be used to quantify errors in the load estimates.


## 2.1. Markov Parameter Representation

Any continuous system may be spatially discretized as ann degrees of freedom (DOF) system, whereby its equilibrium equation of motion can be written in matrix form as
(1)Mx¨t+Cx˙t+Kxt=ft.
For simple systems, M, C, and K 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
(2)u˙t=Acut+Bcft,
where Ac is the system matrix and Bc is the input matrix for continuous case given by
(3)Ac=0I-M-1K-M-1C,Bc=0M-1.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:(4)uti+1=Aduti+Bdfti,
where ti is the subscript over the discretized time and Ad and Bd are related to their continuous case counterparts as
(5)Ad=e[Ac]Δt,Bd=Ac-1Ad-IBc.
Equation (4) is the discrete linear time-invariant state-space representation of the structural system. The corresponding system output is given by
(6)yti=Cduti+Ddfti.
Given that unit impulse load is applied to the system at time t=0, that is, δ0=1, and assuming zero initial conditions, the impulse response at various time points is given by
(7)H0=Cdu0+Ddδ0=Ddu1=Adu0+Bdδ0=BdH1=Cdu1+Ddδ1=CdBdu2=Adu1+Bdδ1=AdBdH2=Cdu2+Ddδ2=CdAdBdu3=Adu2+Bdδ2=Ad2BdH3=Cdu3+Ddδ3=CdAd2Bd⋮Hi=CdAdi-1Bd.
Given the impulse response, the response at any time is given by the convolution of the input force with the impulse response as
(8)yti=∑i=0tiHifti-i,
where the matrices Hi are called forward Markov parameters [20] and summarized as
(9)H0=DdHi=CdAdi-1Bd,i=1,2,3….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 forcesf, given the forward Markov parameters and system responses.From (8), the system response at ti=0 is given by
(10)y0=H0f0.
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 ti=0 is given as
(11)f0=H0TH0-1H0Ty0.
Similarly, at ti=1, the system response is given by
(12)y1=H0f1+H1f0
which can again be solved for the loads as
(13)f1=H0TH0-1H0Ty1-H1f0.
By induction, the loads at timeti can be verified to be
(14)fti=H0TH0-1H0Tyti-∑i=1tiHtifti-i.Next, define a discrete inverse system similar to (4) and (6), where the roles of input forces f and system response y are reversed, as
(15)uti+1=AdIuti+BdIyti,fti=CdIuti+DdIyti,
where the inverse matrices AdI, BdI, CdI, and DdI are related to the corresponding forward matrices by
(16)AdI=Ad-BdDdTDd-1DdTCd,BdI=BdDdTDd-1DdT,CdI=-DdTDd-1DdTCd,DdI=DdTDd-1DdT.
Similar to (8), given the system response, the input forces at any time are given by the convolution relation as
(17)fti=∑i=0tihiyti-i,
where the matrices hi 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 hi are related to the forward Markov parameters Hi by a linear predictive equation given by
(18)h0=H0TH0-1H0T,hk=-h0∑i=1kHihk-i.There exist cases whereH0 is zero matrix. In certain other nonminimum phase structural systems, H0 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]
(19)yti=Cduti+CdBdfti.
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 firstz Markov parameters in (7) are zero, that is, Dd=CdBd=CdAdBd=⋯=CdAdz-2Bd=0, the noncausalz-lead inverse system can be constructed as
(20)hz=HzTHz-1HzT,hk+z=-hz∑i=1kHi+zhk-i+z
and the corresponding force estimates are given by
(21)fti=∑i=0tihi+zyti-i+z.
It must be noted in (21) that the input force estimates f at current timeti are dependent on the structural response y at future times ti-i+z.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 h with the noisy structural response y 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 loadsf can be estimated precisely from measured structural response y 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 zth inverse Markov parameter hz is computed by the least-squares inversion of the zth Markov parameter Hz. All the other higher order inverse Markov parameters depend on hz. Therefore, the accuracy in the computation of hz from Hz dictates the accuracy of all the other higher order inverse Markov parameters. Each row in Hz 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 Hz; the subset so obtained is called the candidate set Hz_cs.To obtain a meaningful solution to (20), an optimal subset Hz_opt of the candidate set Hz_cs needs to be identified, which in turn will lead to an increased accuracy in the estimation of f. As more sensors are used, the additional information helps to obtain a more precise estimate of f, 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 nf, then in order to minimize the error in f estimates, the number of sensorsa must satisfy the criterion a≥nf. Each sensor location, referred to as a candidate point, provides a potential row for inclusion in Hz_opt. Hz_opt needs to be such a subset of Hz_cs that provides the most precise estimates of fti. An approach to extract Hz_opt from Hz_cs is described next.


## 2.4.D-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 fti estimates is given by [22]
(22)varfti=σ2HzTHz-1.
The matrix HzTHz-1 is known as the sensitivity of Hz. For a given variance σ2 in response measurements, minimization of the sensitivity of Hz leads to an increased precision in fti estimates. The sensitivity of Hz 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 Hz leads to the minimization of variation in the fti estimates.For a given number of sensorsa, the candidate set Hz_cs is searched to determinea sensor locations that provide the least variance in the load estimates. One of the criteria followed involves the maximization of the determinant of HzTHz,HzTHz. Designs that maximize HzTHz are called D-optimal designs [23–25], whereD denotes determinant. In order to construct ana-point D-optimal design, thea sensor locations that maximize HzTHz must be selected from the candidate set Hz_cs. To select thea-point D-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 setHz_cs and the number of sensorsa, the first step is to randomly selecta distinct candidate points from the candidate set to initialize Hz. Out of the remaining candidate set, a candidate point is then selected and the corresponding row is augmented to Hz to form Hz+ such that HzT+Hz+ is maximum. Next, out of the a+1 rows in Hz+, a row is deleted to arrive at Hz- such that HzT-Hz- is maximum. This process of adding and deleting rows continues until there is no further improvement in the value of HzTHz. The final Hz so obtained is the D-optimal design Hz_opt 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 obtainM=HzTHz and then calculate the determinant M. This approach of determinant calculation is computationally very expensive. An alternate formula for computing the determinant M+=HzT+Hz+ from that of M when a row yT is augmented to Hz is [25]
(23)M+=M1+yTM-1y,
where + denotes addition. The + is replaced by subtraction for the case when a row yT is deleted from Hz+. In order to be able to use (23), M-1 can be maintained and updated as the row yT is augmented to Hz by
(24)M+-1=M-1-M-1yM-1yT1+yTM-1y,
where - denotes subtraction and is replaced by addition for the case when a row yT is deleted from Hz+.OnceHz_opt 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 f. 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(25)ε%rms=∑fapp-frec2∑fapp2×100%.
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 Figure1. Masses m1 and m15 are connected to fixed boundary. A periodic forcing function is applied to mass m7. The task is to determine the optimum accelerometer locations and reconstruct the input force based on the acceleration time response at those locations.Figure 1
15-Dof spring-mass-damper system.A total of 2 accelerometers were used to measure accelerations of 2 masses. TheM, C, and K 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 setHz_cs. Hz_cs was subjected to the D-optimal design algorithm described in Section 2.4 to obtain Hz_opt whereby the optimum accelerometer locations were found to be at masses m6 and m8. 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 y computed at the optimum accelerometer locations, the input force f 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.Figure 2
Applied and recovered loads at mass 7 with optimum accelerometer placements.


## 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 Figure3. 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.Figure 3
Finite element model of girder with applied load and optimum accelerometer locations.TheM and K matrices were obtained using finite element method. ANSYS provides data for M and K 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 Hz_cs. Hz_cs was subjected to the D-optimal design algorithm described in Section 2.4 to obtain Hz_opt. 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 Hz_cs. 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 y computed at the optimum accelerometer locations, the input force f 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%.Figure 4
Applied and recovered loads using optimum accelerometer placements, no acceleration errors.In order to simulate a more realistic scenario where accelerations are measured experimentally, the acceleration datay 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%.Figure 5
Applied and recovered loads using optimum accelerometer placements, with acceleration errors.Next, a different nature of force was considered where a half-sine wave loading was applied. Both exact acceleration datay 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.Figure 6
Applied and recovered loads using optimum accelerometer placements, no acceleration errors.Figure 7
Applied and recovered loads using optimum accelerometer placements, with acceleration errors.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 ofD-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.



---

*Source: 605695-2015-02-18.xml*