Use of response surface metamodels for identification of stiffness and damping coefficients in a simple dynamic system

Metamodels have been used with success in many areas of engineering for decades but only recently in the field of structural dynamics. A metamodel is a fast running surrogate that is typically used to aid an analyst or test engineer in the fast and efficient exploration of the design space. Response surface metamodels are used in this work to perform parameter identification of a simple five degree of freedom system, motivated by their low training requirements and ease of use. In structural dynamics applications, response surface metamodels have been utilized in a forward sense, for activities such as sensitivity analysis or uncertainty quantification. In this study a polynomial response surface model is developed, relating system parameters to measurable output features. Once this relationship is established, the response surface is used in an inverse sense to identify system parameters from measured output features. A design of experiments is utilized to choose points, representing a fraction of the full design space of interest, for fitting the response surface metamodel. Two parameters commonly used to characterize damage in a structural system, stiffness and damping, are identified. First changes are identified and located with success in a linear 5DOF system. Then parameter identification is attempted with a nonlinear 5DOF system and limited success is achieved. This work will demonstrate that use of response surface metamodels in an inverse sense shows promise for use in system parameter identification for both linear and weakly nonlinear systems and that the method has potential for use in damage identification applications.


Introduction
The concept of a "metamodel" is introduced as a fast running surrogate model.Usually it is a model that will run in minutes on a single processor desktop PC.It could be a traditional reduced order model (as obtained from Guyan reduction or Craig-Bampton methods), a neural network, or a statistically derived model.The focus of this research will be directed toward statistically derived metamodels, or response surface metamodels, for use in system identification applications.Response surface methods may be employed with low computational effort and have the potential to be applied to both linear and nonlinear problems.They also have the advantage of relating input parameters to output features in functional (often polynomial) form, and hence sensitivities are readily identifiable.
Because of the large scale of many damage identification projects, engineers rely heavily on simulation because it is often not feasible to conduct many experiments to explore all damage scenarios.(Imagine an aerospace firm  Fig. 1.Conceptual drawing of the spring-mass system simulated in MatLab(tm) [10].
fabricating and then destroying 500 airplanes to test all possible damage scenarios!)Additionally, with the advent of programs such as the Department of Energy's ASC (Accelerated Strategic Computing Initiative) at National Laboratories across the country, simulations have become increasingly complex, which means they often take tremendous amounts of computing time to run.One example is a simulation of the response of a mechanical joint subject to transient excitation which took 504 processors three hours to run three milliseconds of simulation time on Los Alamos National Laboratory's parallel super computer, Blue Mountain [1,2].In many commercial settings, computing power is simply not available to conduct such complex simulations.While having unlimited experiments and/or simulations would be ideal for exploration of a dynamic phenomenon, it is an unrealistic goal due to budget and time limitations.
The concept of a design space is introduced as the set of all possible experiments or simulations that interest the analyst [3].It is all input variables set at all possible levels and associated dependent features of interest.The design space may be n-dimensional in size.Because the total design space is often prohibitively large, methods have been developed in the literature to efficiently explore it.Through the use of response surface metamodels (and other surrogate modeling methods), analysts and experimenters may be able to efficiently explore the design space to better determine areas of interest and quantify the variability of predictions.
When response surface metamodels are the reduced order model of choice, design of experiments methods (often called response surface methods) are often employed.Already popular in the chemical and industrial engineering communities, design of experiments is a statistical method used to "intelligently" determine which simulation or physical experiments should be run when resources are scarce.The method of design of experiments relies on analysis of variance, or ANOVA, to choose a few points out of the full factorial set that efficiently provide information about the full response space.Metamodels may be fit to the data points using standard multiple regression methods resulting in a polynomial model that relates input parameters to output features [3].While these models are empirical in nature, they rely on the expertise of the experimenter for assignment of model input parameters and choice of appropriate output feature(s).
It is the empirical nature of response surface models that may make them well suited to nonlinear dynamics simulation and experiment, because (depending on the number of training points used) higher order models may be used to relate input variables to response features.In addition to features derived from the frequency domain (which may be difficult to derive in nonlinear settings), response features may be derived from measured or simulated responses in the time domain [4].In fact, the number and types of response features used are limited only by the ingenuity of the experimenter.
To date, structural dynamics researchers have utilized response surface metamodels in structural dynamics for sensitivity analyses [5], uncertainty quantification [6], and model updating [7].However, use of low order models will have an increasingly important part to play in applications such as structural health monitoring and damage  prognosis, where damage identification is considered a fundamental first step [8].Such applications may involve putting many sensors on a structure and processing the data onboard, presumably with some kind of minimal hardware.Because of the low training requirements and simple model form of response surface metamodels, they could be good candidates for use with onboard processing applications, requiring a minimum of computational power.Damage in the structure might be identified by a response surface model and sent to a remote location where an engineer is monitoring structure status.The purpose of this research is to demonstrate on a simple analytical system that response surface metamodels show promise for use in the area of damage identification by demonstrating their ability to identify parameters frequently used as damage indicators.In this paper, response surface metamodels will be fit first for (1) an analytical, linear, five degree of freedom system and (2) an analytical, nonlinear, five degree of freedom system.Only a fraction of the full design space will be utilized to train the models.Once a forward (polynomial) relationship is established and fit accuracy assessed, changed stiffness and damping coefficients will be identified using a standard optimization routine for both the linear and the nonlinear system.Predictive accuracies of the inverse models are assessed.Response surface models are used successfully for parameter identification of stiffness and damping on the simple linear system and with limited success on the simple nonlinear system.

System definition
The system of interest is a simple 5 degree-of-freedom spring-mass system (5DOF hereafter) that has been simulated in MatLab(tm) [10] and shown in Fig. 1.First, stiffness coefficients of the linear system were perturbed (one at a time) with the parameters of mass and damping held constant.Then damping coefficients were changed (mass and stiffness held constant).Once it was established that the location and magnitude of the coefficients could be predicted successfully in a linear system, the procedure was repeated with a nonlinear system.In that system, a nonlinear spring was inserted in place of a linear spring at varying locations, with varying nonlinear stiffness coefficient.This tested whether the method could identify the magnitude and location of a nonlinearity in the system.Next a change in damping coefficient was introduced; but this time the baseline system had a nonlinear spring in place of the linear spring at Location 1.This tested whether the method could identify a changed damping coefficient in the presence of a nonlinearity.The test matrix (Table 1) summarizes the scenarios examined in this research.

Parameter identification of the linear system
Parameter identification of the linear system was first performed for changes to both stiffness and damping coefficients.Stiffness coefficient was chosen because of its familiar use as a damage indicator; most of the time damage to a structural system results in a more flexible (ie reduced stiffness) system.Damping was also chosen for the same reason; it has been recognized that a change in damping always occurs with damage to a system [11].However the difficulty of experimentally or analytically determining damping coefficient has prevented its widespread use as a damage indicator.Use of empirical response surface metamodels may provide an easier way to identify damping parameters.

Identification of stiffness coefficient -Linear 5DOF system
For performance of a response surface metamodel as a system parameter identifier, the relationship between a measureable output feature and the system parameters of interest are of primary importance.In the case of the linear system, relationships between the perturbed stiffness parameter and natural frequencies of the system are physical in nature.Hence the first five natural frequencies were chosen as measurable output features.For the baseline system all spring coefficients had a value of 6000 N/m.Change in the system was defined changing the stiffness coefficient in one of the springs (springs could have a stiffness coefficient of 2000 N/m, 3000 N/m, 4000 N/m, 5000 N/m, or 6000 N/m).The entire design space of interest consists of two input variables, stiffness coefficient and location, each of which could be set at five levels.Hence the full factorial set (all parameters set at all levels of interest) would then consist of 5 levels (2 input parameters) = 25 simulations.
Note that the baseline system was included in the design.The impulse excitation magnitude was 5 N at Location 3. Design of experiments methods were used to choose the input parameter levels that would be used for fitting the response surface metamodel.A three-level design (face centered cubic) was chosen and response surface metamodels were fit to a reduced design space consisting of 9 simulation runs, seen in Table 2. Ideally this set of 9 runs would adequately represent the full factorial set of 25 runs.While 25 runs would not be a computational "hardship" in this simple example, we wanted to demonstrate that the response surface could adequately capture the behavior of the full response space with only a fraction of the full factorial runs.For each of the 9 runs in the design, the first five natural frequencies of the system were retained.Then response surface metamodels were fit using multiple regression techniques.An example of one of these surfaces is shown, in contour format, in Fig. 2.These models relate stiffness parameters (location and magnitude of the coefficient) to the measurable outputs (natural frequencies).One model was fit per natural frequency due to the functional form of polynomial metamodels.A summary of the polynomial models fit is shown in Table 3.It can be seen that the Adjusted R 2 values for all of the models indicates very good fit (see Appendix for calculation of Adjusted R 2 ).

Location Predicted v. Actual, sorted by Stiffness
Once the relationships between the stiffness parameters and each of the five natural frequencies have been established, parameter identification is completed by optimizing stiffness and location parameters using the resulting response surfaces.That is, the question "Knowing the natural frequency values, can the magnitude and location of the reduced stiffness be found?" is answered.Since natural frequencies are not highly correlated (they provide linearly independent information), optimization works well.
The optimization procedure used for this set of models was a Nelder-Meade simplex routine, implemented in Design Expert TM [12].It also incorporates the desirability function of Myers and Montgomery [3] which is as follows: 1 Bolded rows are for the baseline system setup.(1) The powers s and t weight how quickly the function travels to a maximum desirability of 1 at target value B, from the range endpoints of either A or C. Desirability is a criterion placed on the measured output values provided.The optimization routine then changes the value of the input parameters sought until maximum desirability is achieved near the target output feature values provided.The single number desirability criteria is just a multiplication of the desirabilities of each of the output features, in this case the five natural frequencies.Results are then ranked in order of decreasing desirability.In all cases, the standard followed was to choose the set of input parameters that resulted in the highest desirability rating on the output features.
To test how well the optimization procedure worked, the full factorial set of experiments was tested using the parameter identification procedure described above (this includes the 9 design points plus 16 points not included in the design, for a total of 25 points).That is, knowing the five natural frequency values for each of the 25 simulations belonging to the full factorial set (derived from the Matlab(tm) simulation [10]), could the change in stiffness (if any) and the location of the change be identified?Results are shown in Fig. 3 in the form of two plots of actual stiffness or location versus predicted stiffness or location.Ideally these results would plot on a line at 45 • .It can be observed from Fig. 3 that while stiffness predictions follow the correct trend, as the actual values increase, the corresponding predicted values are slightly less.In fact for the baseline condition, the stiffness coefficient was predicted as 5206 N/m, lower than the actual stiffness of 6000 N/m.Error on stiffness predictions was on average about 10%.Location of the changed stiffness typically had an error of one position or less (e.g.Location 2 predicted as Location 3).These results were achieved by training on 9 runs, or 36% of the full factorial design.The ability to predict location and magnitude of stiffness coefficients from frequency data with high accuracy might imply that the method has potential for use in damage identification applications of more complex systems as well.

Damping coefficient identification -Linear 5DOF system
The next scenario examined was the case of changed damping coefficient in the linear 5DOF system.Mass and spring coefficients were constant at 7 kg and 6000 N/m respectively.The system was subject to an impulse load at Location 5 and the resulting response monitored at Location 3.While natural frequencies of the system worked well as output features with stiffness as a damage parameter, they are not as sensitive to changes in damping.Therefore alternate output features were chosen, temporal and spectral moments, which are often used to describe transient dynamic phenomenon (see [4,13] for detailed discussion of the use of moments as extracted features).Generally moment features involve successive integrations of the history of interest multiplied by successive powers of time: where M i is the i th order moment, t is time and f (t) in this case is displacement as a function of time.The first three temporal moments were used.The first temporal (TE, i = 0) represents how much total energy there is in a displacement time history and the second temporal moment (TT, i = 1) is the point in time at which half of the energy has "arrived."The third temporal moment (TD, i = 2) is the RMS duration of the displacement time history.Fig. 6. Results of solving the damage identification problem using optimization.Note all healthy system predictions had a predicted value of Location 4, and are not shown on the location plot (healthy is defined as all dampers set at 0.75 N-s/m).
The first and third spectral moments were also used as features and their description is similar to their analogous time history counterparts.The first spectral moment (SE) represents how much total energy there is in the broadband sense and the third spectral moment (SD) represents the RMS duration of the signal.Figure 4 shows a displacement time history and Fig. 5 a frequency response function for one simulation with data taken at Location 3.While all nine data sets look nearly identical, use of integral moment features allows the analyst to distinguish between them.
The baseline system was defined as all damping coefficients set at 0.75 N-s/m.Damping coefficients were changed one location at a time and could have values of 0.25 N-s/m, 0.50 N-s/m, 0.75 N-s/m (baseline), 1.00 N-s/m, or 1.25 N-s/m.The entire design space of interest consisted of two input variables (location and magnitude of damping Fig. 7. Results from damage identification using response surface metamodels developed for a nonlinear system.(Diamonds are from the design set, triangles are the rest of the points in the full factorial design).coefficient), each of which could be set at 5 different levels for a total of 25 possible simulations in the full factorial design.Note that the baseline system was again included several times in the training set due to the form of the model.Response surface metamodels were fit to a face centered cubic design space (damping and location each set at their "high," "middle," and "low" values) consisting of (3 levels) (2 input parameters) = 9 simulation runs.
For each of these 9 runs, the output features described above were retained (see Table 4).Then response surface metamodels were fit using multiple regression techniques.Polynomials are not shown for brevity.Adjusted R 2 values again demonstrated very good fit.
To test how well the optimization procedure worked, the full factorial set of experiments was tested in the above described parameter identification (optimization) routine (this includes the 9 design points plus 16 points not included in the design, for a total of 25 points).That is, knowing the five output feature values (derived from the Matlab(tm) simulation) for each of the 25 different combinations of damping and location, could the change in damping (if any) and the location of the change be identified?Results are shown below in Figure 6 in the form of two plots of actual damping or location versus predicted damping or location.Ideally the results would plot on a line at 45 • .In Fig. 6 it can be seen that all damping coefficients are predicted well, with the largest error (∼32%) at the lowest damping value.Locations are predicted well; all of the design point locations (1, 3 and 5) are predicted with very high accuracy.Location 2, which was not in the design set, however, is consistently underestimated, causing overlap with the range of predicted values for Location 1.
Most of the errors on damping predictions are under 10% (with the exception of some of the low damping coefficient values).Location error, while high for Location 2, was otherwise generally less than 5%.Higher error was expected at the damping coefficients (0.50, 1.0) and locations (2, 4) that were not included in the training data set and this was observed.Response surface metamodels trained on 36% of the full factorial design space were able to satisfactorily predict damping coefficients and location of the changed values given the simulated output features.Often damping is not known in an absolute sense, but response surface methods could take full advantage of knowledge of damping in a relative sense.The ability to locate and predict magnitude of a damping coefficient in a linear system shows that the method has been used to identify damage-indicating parameters, and that its potential to assess changes in a more complex system should be explored.

Parameter identification of the nonlinear system
Because response surface metamodels were shown to be useful in the context of linear problems, for identification of location and changes in stiffness and damping coefficient identification and location, it was a natural extension to test their ability to detect changes in a nonlinear setting.In this study, perturbations were introduced in the simulation in the form of a nonlinear spring (kx 3 ) that replaced one of the linear springs in the 5DOF spring-mass system.The location and magnitude of the nonlinear stiffness coefficient are then predicted.The second scenario examined was to make the baseline system nonlinear, by replacing the linear spring at Location 1 with a nonlinear spring, and then try to identify a changed damping coefficient and its location in the presence of this nonlinearity.

Introduction of a nonlinear stiffness coefficient -Nonlinear 5DOF system
The baseline was again defined as all linear springs with stiffness coefficient set at 6000 N/m.A change in the system was defined as a nonlinear spring replacing a linear spring at any of the five locations.Degree of change was defined as the value of the nonlinear spring coefficient, which could have the following values: 6e6, 7e6, 8e6, 9e6, 1e7 N/m 1/3 .
The same design was again chosen, a 9 run fraction of the full factorial design space (25 runs), to formulate relationships between controllable inputs and measurable outputs.Because of the presence of nonlinearities in the system, natural frequencies were not chosen, because they were more difficult to identify.Instead maximum displacement (MD) and spectral and temporal moments were chosen as measurable output features [4,13].In this case, the first and second spectral moments (SE, ST) were chosen, along with the first two temporal moments, TE and TT (see previous section for a physical description of moments).The face centered cubic design is again shown in Table 5.Power law transformations of the output features were used to improve fit quality in some cases.Even though the adjusted R 2 statistic for maximum displacement, spectral energy and temporal energy were low (<0.8),their inclusion still yielded better predictions than optimization with just the subset of features with higher adjusted R 2 values.
Results from the parameter identification (optimization) procedure are shown in Fig. 7. Optimizations for the design set and the full factorial set have been plotted separately for both location and stiffness coefficient.Location is predicted well (with one exception) for both the design set and the additional points of the full factorial set.Nonlinear stiffness is predicted reasonably well for the design set (errors of less than 15%).However stiffness predictions for those points outside of the training set are significantly worse (errors over 20% in a few cases).Prediction of the nonlinear stiffness coefficient could be improved with a better choice of output features.However, finding new output features is a difficult task since they should be intuitively linked to the physics of interest, as well as independent of one another so that they do not provide the optimization with redundant information.
For the case of the nonlinear 5DOF system, in which damage is a nonlinear spring with varying stiffness coefficient inserted at one of the five possible locations, response surface metamodels were able to satisfactorily predict the location of the changed stiffness coefficient, but not the magnitude using 36% of the design space.

Damping coefficient identification -Nonlinear 5DOF system
Next, a nonlinear spring (in the form of a kx 3 term, with k = 6e6 N/m1/3 ) replaced the linear spring at Location 1, making the baseline 5DOF spring-mass-damper system nonlinear.Displacements were kept in the same regime as that of the linear spring (see Fig. 8).Damping was again changed at various locations using the same values as for the linear system.The displacement based output features that had been previously proven on the linear system were tried first, however location of the changed damping was not predicted well.In an effort to improve the predictions, velocity based temporal moments, spectral moments and maximum velocity were chosen as features.

Actual
Predicted Ideal Predicted Fig. 11.Results from damage identification using response surface metamodels developed for a nonlinear system.
The face centered cubic design of 9 runs, a fraction of the full factorial design space (25 runs), was again chosen to formulate relationships between system parameters and measurable outputs.Figure 9 shows a sample velocity time history and Fig. 10 shows an FRF plot of a system with reduced damping at Location 3 (0.25 Ns/m).It was found that the nonlinearity at Location 1 had the effect of amplifying changes in damping more than the other locations.In Table 6 it can be seen that the range of output feature values corresponding to the high and low damping values at Location 1 is larger than the output feature ranges for the other locations.
Results from the damage identification (optimization) procedure are shown in Fig. 11.The baseline system had a predicted damping of 0.78 N-s/m.While the two lower values of damping are confounded (same value predicted for both), the higher coefficients are predicted quite well, with errors generally less than 15%.It can be seen that location prediction ranges all span at least one full location value.The response surface method, with the current feature set and the face centered cubic design is capable of predicting a change in damping in a system, but is not able to locate it.Predictions could be improved through a better choice of output features, or potentially an increased number of training sets so that the order of the response surface may be increased.

Conclusions
The method of generating forward relationships with relatively few "training" data sets and the application of the inverse formulation to conduct parameter identification and location with some quantifiable degree of accuracy has been demonstrated on a simple 5 DOF system.Response surface metamodels show promise for application in damage identification (in the form of changed stiffness and damping coefficients) of both linear and nonlinear dynamic structures.They have the advantage of requiring relatively few training data sets as well as the ability to capture parameter changes and locations in some nonlinear problems.Predictive ability of response surface metamodels (as well as other reduced order modeling methods) is highly sensitive to choice of low order feature, or feature extraction, a current topic of research in structural dynamics.Determination of output features remains the single most difficult (and important) task in the process of generating response surface metamodels; sensitivities derived from the response surface in a forward sense should be leveraged when choosing.The order of the response surface models may be increased, but appropriate features must still be chosen and they must represent the physics of interest, as well as provide independent information to an optimization routine.Future work will further investigate the utility of response surfaces for identification of damping parameters, either in an absolute or relative sense.Issues of uncertainty propagation through the inverse process have been addressed in an empirical way in [9], but more research must follow.The methods must be tested on more complex systems than this simple 5DOF analytical example.Finally the efficacy of the response surface method of parameter identification should be compared to other reduced order methods.

Fig. 2 .
Fig.2.One of the five response surfaces generated for the damage identification problem.Model input (damage) parameters are stiffness and location, output in this plot is the second natural frequency.(Plot generated in Stat Ease Design Expert[12]).

Fig. 3 .
Fig.3.Results of solving the damage identification problem using optimization.Healthy system stiffness predictions are not shown, as they were all predicted with a stiffness of 5206 N/m at Location 1.

Fig. 9 .
Fig. 9. Velocity time history for damaged damping at Location 3. Integral features allow the use of data that by eye looks very similar from one damage scenario to the next.Damping at Location 3 = 0.25 N-s/m.

Fig. 10 .
Fig. 10.FRF plot for damaged damping at Location 3. Integral features allow the use of data that by eye looks very similar from one damage scenario to the next.Damping at Location 3 = 0.25 N-s/m.

Table 2
Fractional factorial design used for damage identification of the 5DOF system1 1Bolded rows are baseline system setups.

Table 3
Summary of polynomials generated1

Table 4
Face centered cubic design used for damping damage identification of the 5DOF system1

Table 5
Design used for identification of nonlinear stiffness coefficient

Table 6
Face centered cubic design used for damage identification of the 5 DOF system with nonlinearity and damping damage1