Full text of DOI:10.1007/PL00013285

Struct Multidisc Optim 22, 307–321  Springer-Verlag 2001
Development of simplified models for design and optimization
of automotive structures for crashworthiness
C.H. Kim, A.R. Mijar, J.S. Arora
Abstract Simplified models can be useful for up-front
design of automotive structures for passenger safety dur-
ing crash. Formulations based on the system identifica-
tion approach are presented for development of simplified
models for simulation and design for automotive crash
environment. Numerical crash data available from experi-
ments or simulations are used in the development of such
models. Parametric as well as nonparametric formula-
tions of the problem are investigated. Standard nonlin-
ear programming optimality conditions and methods are
used to solve the resulting nonlinear identification prob-
lem. Simple numerical examples are solved to illustrate
the proposed formulations and methodologies. As a prac-
tical example, the front horn of an automotive structure
is replaced by a single degree of freedom system (SDOF).
Two basis functions that identify the given target data
are studied: Hat functions (piecewise linear) and Cheby-
shev polynomials. Effects of the number of design vari-
ables on the final solution to the problem are investi-
gated. In addition, using the identified SDOF model, re-
design of the front horn to improve its performance is
discussed.
programs have greatly enhanced the state-of-the-art for
simulation of various automotive crash events, such as
frontal crash, side impact, rear crash, offset crash, etc.
These programs use explicit methods for numerical inte-
gration of the equations of motion, which require a very
small time increment to obtain realistic simulations. In
addition, the finite element models are usually very large
having tens of thousands of degrees of freedom. As a re-
sult, the number of calculations for a realistic simulation
is extremely large requiring enormous computing time on
supercomputers.
Design of automotive structures for passenger safety
during crash uses an iterative process where design
changes are made and the structure is re-analyzed for its
response. An evaluation of any design change requires
simulation of the system, which is a nonlinear transient
dynamic problem. The redesign process can be quite te-
dious and time-consuming when full-scale finite element
models are used and the design process is carried out
manually. Such detailed models and simulations are not
useful at the conceptual design stage where quick de-
sign decisions need to be made. Therefore it is desirable
to develop simplified models that can be used at the
conceptual design development stage for quick analysis
and redesign. The design obtained using the simplified
models can then be cascaded to detailed design of the
components.
Key words automotive structures, crashworthiness,
simplified models, optimization, system identification
The purpose of this paper is to explore various
methods and formulations for development of simplified
models for automotive structure and its components.
Considerable data are available from laboratory/field ex-
periments for structural components and the system.
These data can be exploited in the development of the
simplified models. To explore various methods and formu-
lations, development of single degree of freedom models
is considered first. The methods to be explored are gen-
erally classified as system identification techniques. The
basic idea of these approaches is to determine the sys-
tem properties, given the input to the system and output
from it. For linear systems, these techniques have been
developed and used for many years. However, very little
has been done for nonlinear systems, especially under dy-
namic loading environment (Hollowell 1986). Nonlinear
systems are highly complex, requiring systematic ap-
1
Introduction
Numerical methods for simulation of automotive crash
events have been developed for the last about twenty
years. Many finite element programs, such as DYNA3D
and several of its variants, have become available. These
Received March 2, 2000
Revised manuscript received November 14, 2000
C.H. Kim, A.R. Mijar, J.S. Arora
Optimal Design Laboratory, College of Engineering, The Uni-
versity of Iowa, Iowa City, Iowa 52242 USA
e-mail: Jasbir-Arora@uiowa.edu

308
proaches to be developed for their identification. In this
paper, optimization-based approaches are investigated
for identification of nonlinear dissipative dynamic sys-
tems. Different formulations of the problem are proposed
and evaluated.
the mixed method, which combined the hybrid and the
analytical methods. That approach was demonstrated
on two example problems where some structural com-
ponents were modelled as nonlinear spring elements and
others were modelled using the finite elements. The mixed
method was found to give adequate simulation capability
for design of automotive structures for crash environ-
ment. Ni and Song (1986) also formulated and solved
an optimization problem using the mixed simulation
method. Some of the beam components were modelled
using finite elements. Mass of these elements was taken
as the cost function that was minimized. There were two
design variables: wall thickness of a beam component.
Constraints were imposed on Vehicle Crash Severity In-
dex (VCSI), Windshield Residual Crush distance, and the
A-pillar deformation, which was represented by the front
tunnel deformation (this indirectly assesses the steering
column displacement regulated by Federal Motor Vehicle
Safety Standards).
Mahmood et al. (1993) have described in detail a pro-
cedure for rapid simulation and design of the frame of
an automotive structure. They developed a simplified
program, called V-CRUSH, for rapid simulation of the
structure. The program used special collapsible 3-D thin-
wall beam elements and was used to design full front-end
frame for a light truck. The frame was divided into sev-
eral substructures that were designed and tested. Experi-
ments were also performed on the structures. Correlation
between the experimental and simulation results was very
good.
Huang et al. (1995) described Ford’s Energy Manage-
ment System that used CRUSH (Crash Reconstruction
Using Static History) lumped mass modelling capabil-
ity. In that system, the energy absorbing (EA) struc-
tural components were represented by nonlinear springs.
Force-displacement characteristics of the EAs were ob-
tained through static crush tests. Those were input di-
rectly to the program. Dynamic environment of the crash
event was treated by the velocity sensitivity factors (dy-
namic amplification factors). Using the system, barrier
loads and passenger compartment loads were calculated
and compared to the test results in a frontal crash. Cheva
et al. (1996) have also used lumped mass models for
frontal/offset crash studies. They also used design of ex-
periments approach to perform parameter optimization
for crash environment.
Yamazaki and Han (1998) studied crushing energy ab-
sorption of circular and square tubes. Four node shell and
solid finite elements in DYNA3D program were used to
model the tubes. The tubes had a rigid mass attached
to one end. They were crashed into a rigid wall at some
initial velocity. Several cases were solved with different
thickness and radius (width) of the tubes. Initial imper-
fections were introduced in some of the cases. Depending
on the dimensions, axisymmetric, nonaxisymmetric, or
column buckling modes got generated during crushing.
Some of the results were compared with the available ex-
perimental data. Maximum energy was absorbed for the
In Sect. 2, a brief overview of the literature on the sub-
ject is presented. Section 3 defines a general system iden-
tification problem. Section 4 presents parametric identi-
fication of a single degree of freedom system. Section 5
presents nonparametric identification of the single degree
of freedom system. A quadratic programming formula-
tion for the single degree of freedom system is discussed in
Sect. 6. Section 7 considers identification of the front horn
as a single degree of freedom system. Data from simu-
lation with DYNA3D are used in the identification pro-
cess.Hat functions and Chebyshev polynomials are used
to analytically represent the resisting force in the front
horn. Section 8 considers a redesign of the front horn to
improve its energy absorption capability. Finally, Sect. 9
contains discussion and concluding remarks.
2
Overview of literature
2.1
Crashworthiness analysis
Lumped mass models have been used since early 1970’s
for analysis and design of automotive structures for safety
during crash (Kamal 1970; Kamal and Wolf Jr. 1982). In
these models, major nonstructural components are rep-
resented by lumped masses and major deformable struc-
tural components are modelled as nonlinear spring elem-
ents, typically represented with the force-displacement
data obtained from tubes. Lumped parameter models
have been employed for dynamic simulations and oc-
cupant analysis (Bennett et al. (1991). Car to barrier,
car to car, side impact, off-set impact and other crash
events have been investigated. More recently, finite elem-
ent models and a combination of them with lumped mass
models have been utilized for crash simulation (Hollowell
1986).
Ni and Song (1986) described three methods for simu-
lation of automotive structures for crash environment.
The first method was called the hybrid method, which
used a lumped mass model for the structure. The struc-
tural components were represented by nonlinear spring
elements whose force-displacement characteristics were
obtained by static crush tests in the laboratory. Dynamic
correction or amplification factors were used to convert
these curves for the dynamic environment of the crash.
Applications of the method were limited to mostly one-
dimensional problems. The second method was called the
analytical method, which was based on either the limit
analysis of the automotive structure modelled as a space
frame, or the finite element analysis of the structure using
beam and shell elements. The third method was called

309
axisymmetric crush mode. An approximate response sur-
face as a function of the radius (width) and thickness of
the tube was generated using the “design of experiments”
approach. The tube was optimized to maximize the en-
ergy absorption during impact with side constraints on
the dimensions. It was found that the energy absorption
increased with a decrease of the radius (width) and an in-
crease of the thickness. However, if the radius was reduced
too much, the tube could go into a column buckling mode,
which had a dramatic decrease in the energy absorption
capacity. It was also found that the mean impact force in-
creased linearly with the radius (width) multiplied by the
thickness squared. It was also observed that the energy
absorption capacity of the cylindrical tube was slightly
better than that for the square tube for the same mass of
the tubes.
Lust (1992) considered optimization of automotive
structures with constraints based on both linear elas-
tic and crashworthiness loading conditions. A two phase
crashworthiness analysis approach was used. In phase
one, the structural components were analyzed to develop
their force deformation characteristics. These force-de-
formation characteristics were assumed to be scalable
with respect to the wall thickness of the structural elem-
ents, which were treated as design variables in the op-
timization process. A nonlinear approximation was de-
veloped for the crashworthiness constraints. In the second
phase of crashworthiness analysis, a lumped mass model
was used to determine response of the system. Stress, dis-
placement, and frequency constraints were imposed for
linear elastic response. For crashworthiness, constraints
were imposed on VCSI (Vehicle Crash Severity Index)
and total crush distance of the front end. It was found
that a lighter design was obtained with the consideration
of the simultaneous requirement of performance under
linearly elastic condition and the crash condition.
Yang et al. (1994) presented a feasibility study of
using numerical optimization methods to design struc-
tural components for crash. The presented procedure
required several software, which included parametric
modelling (Pro/ENGINEER), automatic mesh gener-
ation (PDA PATRAN3), nonlinear finite element analysis
(RADIOSS), and optimization programs. Both single and
multiple objective formulations were used for numeri-
cal optimization, which resulted in better designs. It was
found that crash optimization was feasible but costly and
that finite element mesh quality was essential for suc-
cessful crash analysis and optimization. To demonstrate
this work a simplified front horn was used that is also
discussed later in this paper.
The nonlinear springs represented by static test data
with some dynamic factors in the forgoing lumped mass
models are no longer accurate for many current configura-
tions. More accurate descriptions of the structural char-
acteristics in a vehicle are required for acceptable predic-
tions. In addition, data observed in a static crush test are
generally different from those observed in a crash event.
Therefore, it is desired to drive the structural character-
istics from the data given during a structure crash test
(Hollowell 1986). This can be posed as a system iden-
tification problem for the nonlinear dynamic structure.
Therefore in the next two subsections, a survey of lit-
erature for system identification problems for linear or
nonlinear structures is presented. Two types of formula-
tions are discussed: parametric system identification and
nonparametric system identification.
2.2
Parametric system identification
System identification methods can be classified on the ba-
sis of their search space:
– parametric methods that search in parameter spaces,
and
– nonparametric methods that search in function spaces.
Basically, parametric methods seek to determine value
of the parameters for an assumed model of the system
to be identified. Schwibinger and Nordmann (1988) pre-
sented a procedure to identify a reduced order model
for torsional vibration of large steam turbine generators.
A 250 DOF finite element model was successfully reduced
to a 13 degree of freedom model. The basic idea was to
minimize the residuals between the eigenvalues of the fi-
nite element model and the reduced order model for the
lowest few modes of vibrations. A weighted least squares
approach was used, and an iterative numerical algorithm
was used to identify stiffness parameters for the reduced
model.
Lingener and Doege (1988) used the system identifi-
cation approach in the time domain to identify stiffness
and damping properties for the supports of a nuclear re-
actor. The vessel was modelled as a rigid body supported
by a six DOF spring-damper system with three transla-
tional and three rotational degrees of freedom. Experi-
mental data about the support system were collected dur-
ing construction of the reactor. Two objective functions
were tried: weighted sum of the squared errors, and an ab-
solute sum of the errors. The second objective function
performed better, and so results with only that function
were presented.
Udwadia and Shah (1976) presented an approach to
identify structures using records obtained during strong
earthquake ground motion. The history-matching prob-
lem was expressed as a minimization of the error be-
tween the measured response and the model response.
The building was modelled as a shear beam, and stiffness
as a function of the height. Any constraints on the prob-
lem were handled using penalty parameters. The steep-
est descent and conjugate gradient methods were used to
solve the unconstrained problem. The sensitivity analy-
sis of the objective functional was performed using the
adjoint variable method, which was judged to be quite
efficient compared to the direct differentiation approach.
The solution was found to be highly dependent on the

310
starting values for the unknown variables (the stiffness).
It was noted that the nonuniqueness of the solution could
result in stiffness that was 20% to 30% different from the
actual stiffness.
When there is strong coupling, the first approximation
worked better. Four example problems were presented:
simple damped linear system, Duffing Oscillator, Van der
Pol Oscillator, and Hysteretic Oscillator. For all exam-
ples, the system was successfully identified using a known
solution, and then tested using another input function.
The approach was also used for multi-degree of freedom
systems where Chebyschev polynomials were replaced by
arbitrary orthogonal functions (Udwadia and Shah 1976).
A method to identify viscoelastic materials characterized
by Volterra integral equations was presented by Diste-
fano and Todeschini (1973). An optimization based ap-
proach was used by Jao et al. (1991) to identify plastic,
viscoelastic and viscoplastic materials modeled using the
endochronic constitutive theory. The problem had several
nonlinear constraints. These material identification prob-
lems are highly nonlinear.
A detailed overview of the literature presented by
Hollowell (1986) revealed that many techniques had been
investigated for identification of linear systems. This was
reasonable at that time since most of the systems were
designed to operate in the linear regime. The review
was divided into two parts: frequency domain approach,
and time domain approach. Much of the literature was
devoted to the frequency domain approach because fre-
quency related data could be measured more accurately.
Sophisticated test equipment and procedures were avail-
able for measuring natural frequencies and mode shapes.
The identification approach determines the mass, stiff-
ness and/or damping matrices for the system. The time
domain data for impulsive loads, resonant testing or
free decay response could also be utilized to identify the
modal parameters. The methods include the least squares
complex exponential method, the polyreference method,
and the Ibrahim time domain method (Leuridan et al.
1985). To identify damping and stiffness of a building sub-
jected to earthquake excitation, a time domain recursive
least square technique that required no matrix inversion
was used by Caravani et al. (1977). The proposed algo-
rithm provided means to account for both the model
uncertainty and the investigators’ confidence in the initial
guess of the parameters.
Hollowell (1986) attempted an identification of lumped
mass models in the time domain by modelling the struc-
tural components as nonlinear springs. In the approach
developed, the stiffness and damping characteristics of
an element were assumed to be separable. In addition,
they were assumed to be piecewise linear over certain
time intervals. The coefficients of the piecewise linear
stiffness and damping characteristics of an element were
treated as the parameters to be identified. Linear and
quadratic programming techniques were used to solve for
these unknown parameters. Three error norms were tried
as cost functions of the optimization problem: L1, L2 and
L∞. In the most solution cases, the error function based
on the L2 norm gave better results. The methodology
was applied to a three DOF model having six nonlinear
springs and two damping elements. Five time increments
were used in the identification process. Constraints were
imposed on the maximum and minimum values of the
damping force. An adaptive approach was developed to
allow for variation in the number of time steps during
each increment of the identification time interval.
2.3
Nonparametric system identification
Nonparametric methods produce the best functional rep-
resentation for structural elements of the system without
a priori assumptions about the deformable elements of
the model. Therefore, these methods make it possible to
eliminate the restriction of forcing the system character-
istics to fit an assumed form such that the system can
not be changed (Masri and Caughey 1979). Automotive
structures and components must be designed to oper-
ate in the linear elastic range during normal operations;
however, they must deform inelastically during crash to
dissipate a large amount of energy. Therefore, nonlin-
ear systems must be treated. Masri and Caughey (1979)
approximated the unknown restoring force for a single
degree of freedom system in terms of Chebyschev polyno-
mials. Then regression analysis was used to identify the
coefficient of the polynomials. Orthogonality property of
these polynomials gave closed form expressions for the
unknown coefficients. The restoring force was assumed to
depend on the deformation (displacements) as well as de-
formation rate (velocities). Two approximations for the
restoring force were tried:
Masri et al. (1993) introduced a procedure based on
the use of artificial neural networks for the identification
of nonlinear dynamic systems. The proposed method was
applied to the damped Duffing oscillator under determin-
istic excitation. The representation of the restoring force
was given by the neural network topology with its weights
instead of Chebyshev polynomials. Therefore identifica-
tion was to find optimum weights to reduce least-squares
error of the restoring force. It was shown that employing
a three-layer neural net with inputs, one output, and 15
and 10 nodes in the first and second hidden layers, respec-
tively, was adequate to characterize the internal force in
the damped Duffing oscillator. The comparison between
neural networks and another nonparametric identifica-
tion method in a previous paper (Masri and Caughey
1979) was discussed.
– Chebyshev polynomials in terms of two variables (dis-
placement and velocity), and
– sum of two terms, one depending on displacement and
the other dependent on velocity.
The second approximation worked well when there was
weak coupling between the displacement and velocity.

311
Argoul and Jezequel (1989) proposed a new inter-
polation procedure which improves identification of the
nonlinear part of lumped parametric systems in areas of
the state space fields where experimental data are in-
sufficient. A modal representation was used and was ap-
proximated by two-dimensional Chebyshev polynomials
while in Masri and Caughey (1979), Masri et al. (1993)
the restoring force was approximated by Chebyshev poly-
nomials. Two examples were presented: Van der Pol os-
cillator and a two degree of freedom model presenting
a Duffing type nonlinear oscillator. The proposed method
requires information on the identification of the linear
behavior of the structure – its pertinent mode shapes –
which allows change to the modal representation and its
dynamic nonlinear response.
Ni et al. (1999) presented a nonparametric identifi-
cation method for nonlinear hysteretic systems. Making
use of the Duhem hysteresis operator, the multi-valued
relationship of hysteretic restoring force in terms of the
displacement and velocity of the phase plane, is mapped
onto two single valued surfaces in an appropriate sub-
space in terms of the state variables–displacement and
hysteretic restoring force. Fitting the surfaces with the
generalized orthogonal or nonorthogonal polynomials in
terms of the displacement and the hysteretic restoring
force identified the functions describing the surfaces. As
an example, a wire-cable isolation system was analyzed
with Taylor series expansion.
s
appr(b, t) be the analytically generated response data,
where b is a vector representing system properties, such
as stiffness, damping and mass. The vector b could also
contain coefficients of the expansion of the system prop-
erties in terms of some known functions; e.g., in non-
parametric system identification, use of sine and cosine
functions, Chebyshev polynomials, splines, etc. Now the
error function between the given data and analytical rep-
resentation can be defined as
e(b, t) = sgiven(t)−sappr(b, t) .
(1)
The system identification problem is to reduce this error
and at the same time satisfy all the constraints on per-
formance of the system. This can be stated as a problem
to minimize a scalar measure of the error function over
the given time interval, as
min E(b) ,
(2)
b
where E(b) can be defined using L1, L2, or L∞ norms of
the error function. Using L1 norm of the error function,
E(b) is defined as
T
ꢀ
E(b) = |e(b, t)| dt ,
(3)
0
where [0, T] is the time interval of interest. This will be
called L1 cost function. Using L2 norm of the error func-
tion, E(b) is defined as
3
Formulation of the system identification problem
T
ꢀ
E(b) = [e(b, t)]2 dt .
(4)
The basic system identification problem is to determine
properties of the system using some known data. The
data could be available from laboratory or field experi-
ments, or from computer simulations. For dynamic sys-
tem identification, the system parameters to be identified
could be mass, stiffness, and damping properties. If the
system properties are represented by parameters, such
as stiffness parameters, damping parameters, mass, etc.,
then the problem is called parametric identification; oth-
erwise it is called nonparametric system identification. In
either case, the system identification problem is defined as
follows.
0
This will be called L2 cost function. Using L∞ norm of the
error function, E(b) is defined as
E(b) = max |e(b, t)| .
(5)
t∈[0,T]
This will be called L∞ cost function. The problem de-
fined in (4) is usually called the least squares problem.
A preliminary investigation for the error representation
showed that the L2 cost functional worked better than
those with L1 and L∞ cost functions. It is shown later
that the L2 cost function has a positive definite Hessian. If
the constraints are linear, then a global minimum point is
obtained for the problem. Therefore, only the formulation
and results with the L2 cost function are presented later
in the paper.
Problem 1. Minimize a measure of the error between
the known data and the analytical data generated by
a functional representation of the system properties sub-
ject to constraints on the response of the system and/or
response related quantities.
3.2
3.1
Constraints
Definition of an error function
The error between the given data and the analyti-
cally generated data can be written in several different
ways. To define this, let sgiven(t) be the given data and
Some constraints may need to be imposed while solving
the system identification problem. For example, the fol-
lowing constraints may be imposed.

312
– Final displacement must be equal to some specified
displacement.
– Maximum force should not exceed a specified value.
– Energy absorbed must be greater than some specified
value.
Φi(τ) = 2τΦi−1(τ)−Φi−2(τ) , for i ≥ 2 ,
Φ0(τ) = 1 , Φ1(τ) = τ .
(11)
Another popular representation is the piece-wise lin-
ear representation. This can be represented in terms of
the Hat function over a subinterval ti−1(= ti −∆t) and
– Acceleration should not exceed a specified value at
any time.
– HIC (Head Injury Criterion) needs to be bounded.
– VCSI (Vehicle Crash Severity Index) needs to be
bounded.
t
i+1(= ti +∆t) as

0
0 ≤ t ≤ ti−1
ti−1 ≤ t ≤ ti ,
ti ≤ t ≤ ti+1
ti+1 ≤ t ≤ T,
,




 t−t
i−1
Note that some of the constraints are time dependent,
and so they must be treated appropriately in the opti-
mization algorithm.

t −t
i
i−1
Φi(t) =
(12)
t
−t
i+1


,
t
−t
i
i+1



0
3.3
where ∆t is a time increment between control points.
Many other basis functions can be used.
Nonparametric identification
In the remaining sections, some of the foregoing for-
mulations will be implemented and evaluated for identifi-
cation of nonlinear dynamic systems.
For nonparametric system identification, the unknown
quantity must be represented analytically. For this pur-
pose some basis functions can be used as follows:
m
ꢁ
4
s
appr(b, t) =
biΦi(t) ,
(6)
Parametric system identification
i=0
where bi are the coefficients of the expansion, m+1 is the
number of terms, and Φi(t) are the basis functions (time-
dependent).
4.1
Introduction
Several basis functions can be used. For example if one
uses Fourier cosine functions, then the basis functions are
given as
To introduce the system identification problem, a linear
damped single degree of freedom system, shown in Fig. 1,
is considered initially. The system is at rest at time t = 0
and is subjected to an initial velocity v0. The equation of
motion and the initial conditions are
iπt
Φi(t) = cos
.
(7)
T
Mx¨+Cx˙ +Kx(t) = 0 ,
Chebyshev polynomials as basis functions are given as
Φi(τ) = cos(i cos−1 τ) , τ ∈ [−1, 1] .
(8)
for 0 ≤ t ≤ T , x(0) = 0 , x˙(0) = v0 .
(13)
The response for the system can be written in the
closed form for under-damped, critically damped, or
over-damped cases (Hibbeler 1978). The following data
are used to generate the response that is used as the
given information in the system identification process:
M = 31.6 kg, K = 20000 N/m, v0 = −2.0 m/s, T = 0.4 s,
Chebyshev polynomials given in (8) satisfy the weighted
orthogonality property as

0
i = j ,
1


ꢀ
Φi(τ)Φj (τ)
1−τ2
π
2
√
dτ =
i = j = 0 ,
i = j = 0 ,
(9)


−1
π
where the weighting function w(τ) is ( 1− τ2)−1/2. Ac-
cording to the definition of Chebyshev polynomials, the
normalized time parameter τ (−1 ≤ τ ≤ 1)is defined
as
t− T
T
2t−T
2
τ =
=
.
(10)
T
2
From this equation, τ = −1 when t = 0 and τ = 1 when
t = T. In numerical calculations, the following equations
can be used to generate Chebyshev polynomials:
Fig. 1 Single degree of freedom system

313
x
x
f
f
ꢀ
ꢀ
C = 700 N s/m. With these data the system is under-
damped, with undamped natural frequency of 4 Hz
(25.16 radians/s).
D = Fdamping dx = Capprx˙given dx =
0
0
T
ꢀ
Capprx˙given2 dt ,
(18)
4.2
0
System identification problem
where xf is the final displacement at t = T. The data
for the problem is taken as KL = 18000, KU = 22000,
CL = 500, CU = 800. The analytical solution for the iden-
tification problem is b1 = 1, b2 = 0, b3 = 1 and b4 = 0.
The SQP algorithm in IDESIGN software is used
to solve the optimization problem (Arora 1989). Note
that with the stiffness and damping representations
given in (14), the equations of motion become nonlinear.
DYNA3D program is called directly through the USER
routines of IDESIGN to integrate the equation of motion.
Starting from the following initial values for the design
variables:
The system identification problem is to determine the
system parameters, stiffness and damping coefficients, to
match the approximated response to the given analytical
response. Let the approximate stiffness and damping be
represented as
K
appr[xgiven(t)] = [b1 +b2xgiven(t)]Kꢀ ,
C
appr[x˙given(t)] = [b3 +b4x˙given(t)]Cꢀ ,
(14)
where Kꢀ and Cꢀ are given constants, and b1, b2, b3 and
b4 are the unknown parameters to be determined. xgiven
and x˙given are the given displacement and the given vel-
ocity data, respectively. The cost function for the system
identification problem is defined as the integral of the
square of L2 norm of the error function between the ap-
proximated acceleration and the known acceleration of
the mass over the time interval T. This error function de-
pends on the parameters bi in Eqs. (14) as
b1 = 15 , b2 = 10 , b3 = 10 , b4 = 5 ,
with
E = 9661.7 ,
the algorithm converges to almost the exact solution in 13
iterations:
b1 = 1.00189 , b2 = 0 , b3 = 1.00004 , b4 = 0 ,
T
ꢀ
ꢂ
ꢃ
2
E(b) =
a
given(t)−aappr(b, t) dt ,
(15)
with
0
E = 9.0277×10−5
.
where aappr(b, t) is the approximated acceleration, which
is computed by integrating (13) with the nonlinear stiff-
ness, the nonlinear damping coefficient in (14), and the
initial conditions in (13). agiven(t) is the given accelera-
tion (calculated using the analytical solution for single
degree of freedom system for demonstration purposes).
Constraints are imposed on the explicit lower and
upper bounds for stiffness and damping parameters as
follows:
At the starting point some constraints are violated by
1273%; however, at the optimum point, no constraint is
violated. Figure 2 shows the history of the cost function
[E, error in (15)] for the optimization process. Note that
the optimal point is obtained after 13 iterations and as
a stopping criterion, the length of the search direction
is required to be less than 0.01. Figure 3 shows the ex-
act and approximated displacement histories. There is
practically no difference between the analytical and the
K
L ≤ Kappr(ti) ≤ KU
,
C
L ≤ Cappr(ti) ≤ CU
,
i = 1, 2, . . . , n ,
(16)
where the subscript i represents a time grid point except
the initial point. In the actual computation, these con-
straints produce 1600 inequalities, since n is 400.
Also, a constraint is imposed on the dissipated energy
as
D ≤ DI ,
(17)
where DI is the input kinetic energy given as 0.5mv02
(63.2 Nm), and D is the energy dissipated by the damper
given as
Fig. 2 Cost function history

314
5.2
System identification problem
The system identification problem is to determine an an-
alytical representation for the force FS(t) using the given
data for the force or the acceleration of the mass. This
data could be available from experiments or computer
simulations. Let the acceleration data be known. The cost
function for the problem is defined as the minimization of
T
ꢀ
E(b) = [e(b, t)]2 dt ,
(20)
0
Fig. 3 Displacement–time histories for target data and iden-
tified data
where e(b, t) is given as
e(b, t) = agiven(t)−aappr(b, t) .
(21)
approximated histories. Similar matches are obtained for
velocity and acceleration. The active constraints at opti-
mum point are the upper limits for stiffness. Although the
formulation identifies the linear system well, it needs dis-
placement, velocity and acceleration data. This may not
be desirable in some practical applications. It is possible
to formulate and solve the identification problem when
only the acceleration data are known.
Here agiven(t) is the given acceleration data, aappr(b, t) is
the approximate representation of the acceleration, and b
is a vector of unknown parameters treated as design vari-
ables of the system identification problem. To calculate
the approximate acceleration from (19), the force FS(t) is
approximated as
m
ꢁ
FSappr(b, t) =
biΦi(t) ,
(22)
5
i=0
Nonparametric system identification
where bi are the parameters of the expansion, Φi(t) for
0 ≤ t ≤ T are some known shape functions, and m+1 is
the number of terms used to represent the force. Now the
system identification problem is to determine the param-
eters bi to minimize the error function defined in (20).
Constraints may be imposed on the parameters them-
selves, and/or other quantities such as displacements.
The constraint on the force is expressed as
5.1
Introduction
In this section, a nonparametric representation for the re-
sisting force of the structural element is considered. This
representation is more general than the parametric form.
The system is shown in Fig. 4. The equation of motion for
the system is written as
FSL ≤ FSappr(b, t) ≤ FSU
,
(23)
Mx¨+FS(t) = 0 ,
where FSL and FSU are the lower and upper limits on
the force.
for
0 ≤ t ≤ T , x(0) = 0 , x˙(0) = v0 .
(19)
To solve the foregoing optimization problem, an SQP
algorithm available in the program IDESIGN is used.
The program requires evaluation of problem functions
and their gradients. The program DYNA3D is coupled
directly to IDESIGN to evaluate response of the system
for any values of the design variables. The gradients of
the cost function and the constraints are evaluated using
the central finite difference option available in IDESIGN.
The following data for the system are used: M = 31.6 kg,
v0 = −2.0 m/s, T = 0.1 s, FSL = −1600 N, FSU = 200 N,
spring constant k = 20 000 N/m, and the damping coeffi-
cient c = 700 N s/m. These data are used to calculate the
acceleration of the mass that is taken as the given data
in (21). Five points corresponding to the time 0 s, 0.025 s,
0.050 s 0.075 s, and 0.1 s, are used to represent the force;
i.e., m = 4 in (22). In the present example, the Hat func-
tion defined in (12) is used. Therefore the parameter bi
Fig. 4 Nonparametric single degree of freedom system

315
becomes the force at the i-th point. In the numerical so-
lution process, 101 grid points (100 intervals) are used in
the simulation and the integral evaluations. The bounds
on the force shown in (23) are imposed at each time grid
point, and so there are 202 constraints for the problem.
Figure 5 shows estimated acceleration and the target
acceleration curves. It is seen that there is practically
no difference between the two curves. Since the restor-
ing forces are proportional to acceleration from (19), their
histories match quite well also. The displacement and vel-
ocity histories also matched very well. The cost function,
calculated from (20), is almost zero. The values for the 5
design variables are
the base model. Other data are the same except that FSL
is changed to −1800 N.
The cost function is changed to minimize the max-
imum acceleration of the system over the time interval 0
to T,
ꢄ
ꢅ
min max |a(t)|
.
(26)
0≤t≤T
This min-max problem is converted to the standard one
by introducing an artificial design variable amax and an
additional constraint as follows:
min amax
subject to an additional constraint
|a(t)| ≤ amax 0 ≤ t ≤ T .
,
(27)
b0 = −1449.2 , b1 = −1350.5 , b2 = −890.4 ,
b3 = −417.7 , b4 = 50.2 .
,
(28)
At the optimal point, a few of the largest values for the
approximated force are active constraints.
This constraint is also imposed at each time grid point,
just as the force constraints in (23). Thus there are
304 constraints for the problem including the constraint
of (24).
Five terms are used for Hat functions in (22) to
represent the force-time curve as in the previous sec-
tion. All other data are the same as for the original
case. Figures 6–10 show various response histories for
the re-designed system. It is seen that the objective of
Fig. 5 Acceleration–time histories for target data and iden-
tified data
5.3
Redesign of the system
The system identified in the previous section is redesigned
to meet the desired system response. It is required that
change of kinetic energy (D) during the time interval T be
at least 10% more than that of the base system (Db). The
increase in the kinetic energy change yields more energy
absorption for the system. This constraint is expressed as
Fig. 6 Restoring force–time histories for redesign problem
D ≥ 1.1Db ,
(24)
where D and Db are given as
1
1
D = M[v02 −v2(T)] , Db = M[v02 −vb2(T)] .
(25)
2
2
Here vb(T) is the base model velocity at time T (=
0.66983 m/s), and v(T) is the final velocity of the re-
designed system. Note that (24) implies that the final
velocity of the redesigned system is reduced compared to
Fig. 7 Kinetic energy–time histories for redesign problem

316
6
Quadratic programming formulation with L2 cost
function
System identification problem defined in Sect. 5.2 is for-
mulated as a quadratic optimization problem. The cost
function for the optimization is given from (20) and (21)
as
T
ꢀ
ꢂ
ꢃ
2
min E(b) =
a
given(t)−aappr(b, t) dt ,
(29)
0
where agiven(t) is the given acceleration, T is a termina-
tion time for the simulation, aappr(b, t) is the approxi-
mated acceleration obtained from the basis function rep-
resentation of the force element. Substituting the basis
function representation for the force element in (22) into
the equation of motion (19), the approximated accelera-
tion is given as
Fig. 8 Displacement–time histories for redesign problem
m
ꢁ
1
1
a
appr(b, t) = − FSappr(b, t) = −
Φibi =
M
M
i=0
1
−
ΦT b .
(30)
M
Substitute (30) into the cost function in (29) and simplify
it to obtain
1
E(b) = bT Qb+cT b+constant,
(31)
Fig. 9 Velocity–time histories for redesign problem
2
where
T
ꢀ
2
M2
Qij
=
ΦiΦj dt , i, j = 0, 1, 2, . . . , m ,
(32)
(33)
(34)
0
T
ꢀ
2
ci =
a
given(t)Φi(t) dt , i = 0, 1, 2, . . . , m ,
M
0
T
ꢀ
constant = [agiven(t)]2 dt .
0
Equality constraints are imposed on the final displace-
ment and/or the final velocity:
Fig. 10 Acceleration–time histories for redesign problem
x
appr(b, T)
gd(b) =
gv(b) =
−1.0 = aT b+d = 0 ,
−1.0 = sT b+r = 0 ,
(35)
(36)
xf
10% increase in the kinetic energy change is achieved
(D = 61.75 kg m2/s2; Db = 56.11 kg m2/s2). This implies
that the absolute velocity of the mass is reduced from
0.6683 m/s to 0.302577 m/s compared to the base sys-
tem giving more absorption of energy. Note that the
force-time curve for the re-designed system is below
that for the base model. The reason for this behav-
ior is that no constraint is imposed on the displace-
ment of the system, and so the displacement of the re-
designed system has increased from about 0.023 m to
0.082 m. The constraint in (24) is active at the optimal
point.
x˙appr(b, T)
vf
where xf and vf are the given final displacement and
velocity, and xappr(b, T) and x˙appr(b, T) are the approx-
imated final displacement and velocity obtained by inte-
grating the equation of motion. Equations (35) and (36)
can be written as
ꢆ
ꢇ
ꢆ
ꢇ
gd(b)
gv(b)
aT b+d
sT b+r
g(b) =
=
= AT b+d = 0 ,
(37)

317
Fig. 11 Shape of front horn with notches
where A = [a s](m+1)×2 and dT = [d r]. Detailed ex-
pressions for matrix Q, and vectors c, a, s, etc., with Hat
and Chebyshev basis functions are given by Kim et al.
(1999).
The identification problem defined in (31) and (37)
can be solved in a closed form. The Lagrangian function
for this QP problem is written as
Then the simulation data are used to identify the system.
Three objective functions based on L2, L1 and L∞ norms
of the error function introduced in Sect. 3 are investigated
for the identification problem; however, results with only
L2 cost function are presented.
Figure 11 shows overall dimensions of the front horn
used as the sample problem. Figure 12 shows the crash en-
vironment for the system. The data for the problem are
as follows: initial velocity is v0 = 525 in/s (approximately
30 mph), mass of the horn is 0.0114 lbf·s2/in, attached
mass is 3.69 lbf·s2/in (corresponds to half of 1300 kg, the
vehicle weight). The system is travelling from right to
left and crashes into a fixed barrier. To obtain the re-
sponse histories, a finite element model for the horn is
developed and analyzed with DYNA3D program. The fi-
nite element model consists of 1184, 4-node quadrilateral
shell elements.
1
L = bT Qb+cT b+constant+vT (AT b+d) ,
(38)
2
where v is the Lagrange multiplier vector. The optimality
condition gives
Qb+c+Av = 0 .
(39)
In addition, the equality constraints given in (37) must be
satisfied at the optimal point. Substituting (39) into (37),
the Lagrange multiplier vector is given as
v = (AT Q−1A)−1(−AT Q−1c+d) .
(40)
The optimal values of b are computed by substituting
(40) into (39). It can be shown that the Hessian matrix is
positive definite. Since the constraints in (37) are linear,
the QP problem is convex. Thus the solution obtained for
the problem is a global minimum.
The forgoing quadratic formulation with Hat func-
tions and Chebyshev polynomials is applied to the ex-
ample problems for identification and redesign in the next
two sections.
Fig. 12 Crash environment for front horn
7
Numerical examples: Front horn system
The model, shown in Fig. 13, has 7104 degrees of free-
dom. A 40 ms simulation is performed. The following data
are used: material is elastic-plastic strain hardening, yield
stress is 30000 lb/in2, density is 7.1×10−4 (lbf ·s2/in)
/in3, tangent modulus is 60000 lb/in2, Young’s modu-
lus is 3.0 × 107 lb/in2), Poisson’s ratio is 0.3. Figure 14
shows the deformed configuration of the component. At
T = 0.04 s, x(T) = 19.6 in, and v(T) = 456.7 in/s. The
system has not come to rest yet. The kinetic energy of the
system has reduced from 510000 lb·in to 386000 lb·in.
Among many components, the front horn is an important
part of the automotive structure. Design of this compon-
ent can affect performance of the vehicle during a crash.
Therefore it is important to have a good design for this
component. In this section, the system identification
problem for the front horn to reduce it to a single degree
of freedom system as shown in Fig. 4 is addressed. First
a finite element simulation of the front horn is performed.

318
Fig. 13 Finite element model of front horn
Fig. 15 Comparison of target data and identified data using
Hat basis functions
results from IDESIGN and the QP problem match per-
fectly. In the QP problem, 25 design variables are used,
since the first time control point [t0 in (12)] does not start
from zero.
As another basis function to approximate the given
force data Chebyshev polynomials are used. Numerical
results from IDESIGN and QP formulation are given
in Fig. 16. It is seen that Chebyshev polynomials give
Fig. 14 Deformed shape of front horn at t = 0.038 s
All the responses are obtained using the 100 Hz low
pass filter options in LS-TAURUS, the post-processor for
DYNA3D.
Note that in the simulation process, energy absorp-
tion by other components of the vehicle is ignored. Thus
all the input energy is being absorbed due to crushing of
the front horn. This is not a realistic assumption, prac-
tically. However, the generated simulation data are used
in the identification process to develop and demonstrate
the procedure. The error function e(b, t) of (21) is defined
using the force data and an analytical force representa-
tion of (22). m is taken as 25 and the basis functions are
taken the Hat functions and Chebyshev polynomials in
(22). Thus there are 26 design variables. A constraint is
imposed on the final displacement to obtain the desired
response for the identified system, as x(T) = xf , where
x(T) is the final displacement of the identified system and
xf is the specified final displacement, which is given from
the simulation results as 19.6 in.
Fig. 16 Comparison of target data and identified data using
Chebyshev polynomials
The problem is solved using the SQP algorithm in
IDESIGN program Arora (1989). The starting value for
the force is set to 6000 lbs for the entire time interval.
An optimal solution is obtained after 30 iterations. In
addition, the analytical solution of the QP problem is
obtained using a program written in MATLAB with 25
design variables. The target force and two identified force
data set using Hat basis functions, are shown in Fig. 15. It
is seen that the three force curves match quite well. The
Fig. 17 Minimum cost value vs. number of design variables
with Hat functions

319
Fig. 18 Minimum cost value vs. number of design variables
with Chebyshev polynomials
Fig. 20 Redesigned displacement by Hat functions
Fig. 21 Redesigned velocity by Hat functions
Fig. 22 Redesigned acceleration by Hat functions
a little bit smoother curves. However, as the number of
design variables (i.e., the number of terms in the basis
function expansion) is increased, the optimum cost func-
tion value increases with Chebyshev polynomials com-
pared to those with the Hat function. This is seen in
Figs. 17 and 18. This is due to the known characteristics
of the Chebyshev polynomials that give relatively large
errors at the end points. The Hat function representa-
tion does not have this limitation. Thus, it appears that
Hat function representation will work better for practical
applications.
8
Redesign with simplified models
So far system identification has been discussed, which
was to obtain mathematical expressions to represent the
given tabular target data. Once the system is identified,
it is usually required to redesign it to improve its per-
formance. Therefore, a redesign of the horn structure is
considered next. The procedure for redesign of the system
is the same as for the identification problem. It is required
to reduce the final velocity from 456.7 in/s to 430 in/s
keeping the final displacement same as before, i.e., 19.6 in.
The reduced velocity yields smaller kinetic energy at the
final time. Thus more energy absorption would be ob-
tained with crushing of the front horn. It is noted how-
ever that this reduced velocity is chosen arbitrarily for
demonstration of the methodology. Other more mean-
ingful constraints and data will need to be considered
for practical applications. Figure 19 shows that the op-
timum solutions of the redesign problem from IDESIGN
and the QP formulation match each other well. The fig-
ure also shows difference between the new force curve
and the base data (the target data used in the identifi-
Fig. 19 Redesigned restoring force by Hat functions

320
cation problem). The new force curve can now be used
for detailed redesign of the front horn for its improved
performance. The predicted displacement, velocity and
acceleration using Hat functions and QP formulation are
given in Figs. 20–22, with the base data. Figure 21 shows
that the final velocity is reduced from 456.7 in/s to the
desired value (430 in/s).
tural and Multidisciplinary Optimization (WCSMO-3), held
at the University at Buffalo, Amherst, New York, May 17–21,
1999. Support provided by Ford Motor Company for this re-
search under URP with Dr. C.C. Wu as project monitor is also
acknowledged.
References
9
Arora, J.S. 1989: IDESIGN user’s manual. Version 3.5.2.
Iowa City, IA, USA: Optimal Design Laboratory, College of
Engineering, The University of Iowa
Discussion and conclusions
The problem of mathematically representing highly com-
plex nonlinear dissipative dynamic systems with sim-
plified models is addressed. Basically two identification
methods are discussed: parametric identification and
nonparametric identification. In parametric identifica-
tion, the stiffness and the damping constants that are
represented by the linear functions of the displacement
and the velocity of a SDOF system are identified using
the optimization formulation. DYNA3D is used to per-
form integration of the nonlinear equation of motion for
the system. This formulation needs all displacement, vel-
ocity and acceleration data. Secondly, in nonparametric
identification, two topics are discussed, system identifi-
cation and its redesign. Two types of formulations are
introduced:
Argoul, P.; Jezequel, L. 1989: Improvement of a nonpar-
ametric identification procedure used in nonlinear dynamics.
J. Appl. Mech., ASME 56, 697–703
Bennett, J.A.; Lust, R.V.; Wang, J.T. 1991: Optimal de-
sign strategies in crashworthiness and occupant protection.
In: Crashworthiness and occupant protection in transportation
systems, AMD-Vol. 126/BED-Vol. 19, pp. 51–66. New York:
ASME
Caravani, P.; Watson, M.L.; Thompson, W.T. 1977: Recur-
sive least-squares time domain identification of structural par-
ameters. J. Appl. Mech., ASME 44, 135–140
Cheva, W.; Yasuki, T.; Gupta, V.; Mendis, K. 1996: Vehicle
development for frontal/offset crash using lumped parameter
modeling. SAETrans. , Paper No. 960437, pp. 53–58
– a general nonlinear programming formulation, and
– a QP formulation.
Distefano, N.; Todeschini, R. 1973: Modeling, identification
and prediction of a class of nonlinear viscoelastic material (I).
Int. J. Solids Struct. 9, 805–818
Hat function and Chebyshev polynomials are used to ap-
proximate the given data. First, two formulations using
Hat functions and Chebyshev polynomials are demon-
strated for a linear system. Then as a practical example,
a simplified SDOF system for the horn structure of a ve-
hicle is identified and redesigned using the proposed pro-
cedures. Horn structure goes through large displacements
and plastic collapse (buckling). It is shown that a single
degree of freedom model using appropriate shape func-
tions can capture this complex behavior of the struc-
ture quite well. From the study about the effects of the
number of design variables on the optimal solution, it
is observed that the Hat function representation is su-
perior. Once the system is identified, a redesign of the
structure is carried out to improve its energy absorption
capacity.
Hibbeler, R.C. 1978: Engineering mechanics: dynamics. New
York: Macmillan
Hollowell, W.T. 1986: Adaptive time domain, constrained sys-
tem identification of nonlinear structures. Proc. Symp. on Ve-
hicle Crashworthiness Including Impact Biomechanics, AMD-
Vol. 79/BED-Vol. 1, pp. 105–123. ASME
Huang, M.; Chen, R.; Margolin, B. 1995: Use of an advanced
CRUSH model in estimating vehicle impact loading and en-
ergy absorption. In: Crashworthiness and occupant protec-
tion in transportation systems, AMD-Vol. 210/BED-Vol. 30,
pp. 273–285. New York: ASME
Jao, S.Y.; Arora, J.S.; Wu, H.C. 1991: An optimization
approach for material constant determination for the en-
dochronic constitutive model. Comp. Mech. 8, 25–41
Using the identified simplified model, a new force-
displacement curve for the structure can be obtained to
improve its performance. The modified force curve can
then be used for detailed redesign of the structure to ob-
tain desired response characteristics.
The proposed formulations have also been extended to
multi-degree of freedom systems, where several structural
elements need to be identified simultaneously.
Kamal, M.M. 1970: Analysis and simulation of vehicle-to-
barrier impact. SAETrans. 79, 1453–1467, SAE Paper
No. 700414
Kamal, M.M.; Wolf Jr., J.A. (eds.) 1982: Modern automotive
structural analysis. New York: Van Nostrand Reinhold Co.
Kim, C.H.; Mijar, A.R.; Arora, J.S. 1999: Development of
simplified models for design and optimization of automotive
structures for crashworthiness. Technical Report No. ODL-
99.12, Iowa City, IA, USA: Optimal Design Laboratory, Col-
lage of Engineering, The University of Iowa
Acknowledgements This paper is based on a presentation
made by the authors at the 3-rd World Congress of Struc-

321
Leuridan, J.M.; Brown, D.L.; Allemang, R.J. 1985: Time
domain parameter identification methods for linear modal
analysis: A unifying approach. Proc. ASMEDesign En-
gineering Technical Conf. (held in Cincinnati, OH). Paper
No. 85-DET-90
Ni, C.-M.; Song, J.O. 1986: Computer-aided design analysis
methods for vehicle structural crashworthiness. Proc. Symp.
on Vehicle Crashworthiness Including Impact Biomechanics,
AMD-Vol. 79/BED-Vol. 1, pp. 125–139. ASME
Ni, Y.Q.; Ko, J.M.; Wong, C.W. 1999: Nonparametric iden-
tification of nonlinear hysteretic systems. J. Eng. Mech. 125,
206–215
Lingener, A.; Doege, S. 1988: Estimation of dynamic spring
and damping parameters of the supports of a nuclear reac-
tor by means of an adaptive method in time domain. In:
Natke, H.G. (ed.) Application of system identification in en-
gineering. CISM Course and Lecture No. 296, pp. 511–523.
Vienna: Springer
Schwibinger, P.; Nordmann, R. 1988: Optimization of a re-
duced torsional model using a parameter identification pro-
cedure. In: Natke, H.G. (ed.) Application of system identifi-
cation in engineering (CISM Course and Lecture No. 296),
pp. 525–542. Berlin, Heidelberg, New York: Springer
Lust, R. 1992: Structural optimization with crashworthiness
constraints. Struct. Optim. 4, 85–89
Mahmood, H.F.; Wheatley, D.G.; El-Bkaily, M.A.; Bak
kar, J.Y. 1993: On the front end design of automotive vehicle
for crashworthiness. Crashworthiness and occupant protec-
tion in transportation systems, AMD-Vol. 169/BED-Vol. 25,
pp. 311–318. New York: ASME
Udwadia, F.E.; Shah, P.C. 1976: Identification of structures
through records obtained during strong earthquake ground
motion. J. Eng. Industry, ASME, 1347–1362
Yamazaki, K.; Han, J. 1998: Maximization of the crushing en-
ergy absorption of tubes. Struct. Optim. 16, 37–46
Masri, S.F.; Caughey, T.K. 1979: A nonparametric identi-
fication technique for nonlinear dynamic problems. J. Appl.
Mech., ASME 46, 433–447
Yang, R.J.; Tseng, L.; Nagy, L.; Cheng, J. 1994: Feasibility
study of crash optimization. In: Gilmore, B.J.; Hoetzel, D.A.;
Dutta, D.; Eschenauer, H.A. (eds.) Advances in design au-
tomation, DE-Vol. 69-2, pp. 549–556. ASME
Masri, S.F.; Chassiakos, A.G.; Caughey, T.K. 1993: Identifi-
cation of nonlinear dynamic systems using neural networks.
J. Appl. Mech., ASME 60, 123–133