Pergamon
Int. J. Non-Linear Mechamcs, Vol. 31, No. 5, pp. 601 -617, 1996
Copyrght
,?, 1996 Elsevier Science Ltd
Prmted in Great Bntain. All nghts reserved
002&7462/96
$15.00 + 0.00
SOO20-7462(96)00025-X
SIMULATION
Kurtis
Department
OF A CLASS OF NON-NORMAL
PROCESSES
R. Gurley,*
of Civil Engineering
Ahsan
Kareem
and Michael
and Geological Sciences, University
IN 46556, U.S.A.
RANDOM
A. Tognarelli
of Notre Dame, Notre Dame,
Abstract-This
study addresses
the simulation
of a class of non-normal
processes based on
measured samples and sample characteristics
of the system input and output. The class of nonnormal processes considered here concerns environmental
loads, such as wind and wave loads, and
associated structural responses. First, static transformation
techniques are used to perform simulations of the underlying Gaussian time or autocorrelation
sample. An optimization
procedure
is
employed to overcome errors associated with a truncated Hermite polynomial transformation.
This
method is able to produce simulations
which closely match the sample process histogram, power
spectral density, and central moments through fourth order. However, it does not retain the specific
structure of the phase relationship between frequency components,
demonstrated
by the inability to
match higher order spectra. A Volterra series up to second order with analytical kernels is employed
to demonstrate
the bispectral matching made possible with memory models. A neural network
system identification
model is employed for simulation of output when measured system input is
available, and also demonstrates
the ability to match higher order spectral characteristics.
Copyright c~ 1996 Elsevier Science Ltd.
Keywords:
bispectrum
non-normal,
simulation,
random
wind pressure, ocean waves
processes,
neural
networks,
higher
order
statistics,
INTRODUCTION
The complete analysis of dynamic system reliability necessarily includes a statistical analysis of
extreme response. Often, the response of a system under consideration is non-Gaussian
due to
non-normal
input, non-linear system properties, or a combination
of both. The presence of
non-linearities
leads to extreme response statistics that no longer resemble those extreme
models based on Gaussian processes. The importance of the extreme response to system
reliability has prompted much research in the development of techniques to predict these
extreme statistics (e.g. solution strategies for Volterra systems). In order to validate these
extreme prediction models, time domain response simulation is attractive, since the equations of
motion may be integrated directly to include the full non-linearities. The simulation of Gaussian
random processes is well established [l-4].
Progress in the simulation of non-Gaussian
processes has been elusive, but necessary for time domain simulation of system response to
non-Gaussian
input (e.g. large amplitude waves on offshore platforms, and wind pressure
fluctuations on cladding components). This work considers several techniques to simulate
non-Gaussian
stationary random processes concerning wind and wave related processes [S, 61.
The focus in this work is on the transformation
of Gaussian simulations to non-normal
processes, based on information
provided in samples of the desired non-Gaussian
process.
We concentrate
on the class of non-Gaussian
processes typical of localized wind pressures
as well as that associated with the response of non-linear offshore systems to wind and wave
fields. There may exist practical classes of non-Gaussian
processes for which the tools
presented herein are not necessarily appropriate.
STATIC
TRANSFORMATION
METHODS
Probability transformation
Static transforms relating a non-Gaussian
process with its underlying Gaussian process
have been the basis of a variety of non-normal
process simulation techniques. A sample of
*Author
to whom correspondence
should
be addressed.
601
602
K. R. Gurley
et al
static transformation
techniques can be found in refs [7-lo]. The few studies in this context
have looked at simulation based on a target power spectral density and target probability
density function [ll, 121. A summary of several techniques including the use of filtered
Poisson h-correlated processes and z-stable processes is found in a recent book by Grigoriu zyxwvutsrqpo
c131.
An approach
used by Yamazaki and Shinozuka
[12] begins with the simulation
of
a Gaussian process u(t) which is then transformed to the desired non-Gaussian
process y(t)
through the following mapping:
y(t) = F;
‘{aqu)).
(1)
A similar concept utilizing the translation
process has been introduced
by Grigoriu
[7]. Yamazaki
and Shinozuka
use an iterative
procedure
to match
the desired
target spectrum
by updating
the spectrum of the initial Gaussian
process, since the
non-linear
transformation
in eqn (1) also modifies
the spectral
contents.
This
iterative procedure does not guarantee convergence for all classes of non-linear
processes.
For some y(t) there may be no corresponding
Gaussian
form u(t) with a matching
spectrum.
Correlation distortion
The necessity for an iterative procedure may be eliminated if one begins with the target
spectrum or autocorrelation
of the non-Gaussian
process and transforms it to the underlying correlation
of the Gaussian
process. This approach
is referred to as the correlation-distortion
method in stochastic systems literature
[lo, 14,151. For a given static
single-valued
non-linearity
x = y(u), where u is a standard normal Gaussian process, the
desired autocorrelation
of x in terms of y can be expressed as [lo]: zyxwvutsrqponmlkjihgfedcbaZYXW
k=O
where pXXis the normalized autocorrelation
kth Hermite polynomial
given by
Hk(a) = (-
of the non-Gaussian
l)*exp($)$[exp
process, and Hk(U) is the
(-;)I.
An alternative
to the preceding approach is to express x as a function of a polynomial
whose coefficients
are determined
by a minimization
procedure,
e.g. Ammon
[ll].
Alternatively
one may use translational
models involving Hermite moment transformation models described earlier. In this study, we utilize a Hermite model for its convenience and availability
in the literature. A simulation
based on the schematic shown in
Fig. 1 would eliminate the spectral distortion
caused by the non-linear
transformation,
since its inverse is employed
to reverse the distortion.
The simulation
algorithm
is as follows: (i) estimate the auto-correlation
of the mean-removed
normalized
sample
to the autocorrelation
non-Gaussian
process to be simulated
(R,,(z)); (ii) t ransform
of the underlying
Gaussian
process (R,,(z)) by solving for R,,(z) [S] in the following
equation:
R,,(z) = a2CR,,(r) + 2fi:R:,,(r)
Fig. 1. Schematic
of the correlation
+ zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQP
6@Ri,(41,
(4)
distortion
method
Class of non-normal
I”
random
603
processes
loI
200
3w
4w
500
600
700
0 zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
100
2W
300
400
5W
600
7W
600
9x
""0
6W
9x2
measured data
measured data
-
simulation
simulation
Fig. 2. Measured wind pressure signal (top left), a correlation distortion simulation (top right), and
power spectral density and pdf of the measured data and ensemble of 100 simulations.
where
h, =
Y3
4+2Jm’
L
4
=;I
+
w,
-
1
’
Cd=
18
1
J1
+ 2ff: + 6K; ’
and y3 and y4 are the skewness and kurtosis of the fluctuating
process; (iii) simulate
a Gaussian process using the spectrum, G,=(w), associated with R,,(T); (iv) transform this
simulated process u, back to a non-Gaussian
process using
x = a[u + &(u2 - 1) +
h4(u3
-
3u)];
(5)
(v) replace the mean and variance of the original parent process to produce a simulation,
x, of the original non-Gaussian
process x.
Figure 2 compares a measured wind pressure signal with a single realization of a correlation distortion simulation of that signal in the top left and right figures, respectively. Note
the undesirable positive extreme behavior in the simulation that is not seen in the sample.
The bottom figures compare the power spectral density and PDF of the measured data with
that from the ensemble average of 100 correlation
distortion
simulations.
The statistical
moments of the standardized
original signal are compared
with the average moment
statistics of the 100 realizations in Table 1. The higher positive kurtosis in the simulations
can be observed in both the time history and the PDF comparison. Unless otherwise noted,
the comparison
of the statistical properties of the simulation
with those of the sample
process is made using an ensemble of 100 realizations for the sake of expedience. Ensembles
K. R. Gurley
604
Table 1. Statistics
ofmeasured
wind pressuredata
et a/
and ensemble averaged
Standard
deviation
Measured wind data
Ensemble of 100 correlation
distortion simulations
Ensemble of 100 modified direct
transformation
simulations
Coefficient of
skewness
simulated
data
Coefficient
kurtosis
of
1.0
0.9927
- 0.8309
- 0.7960
4.9940
5.6711
0.9960
-0.8120
4.7676
of up to 1000 were used in the initial phases of this work, and showed little difference in the
results. A later example using 2000 realizations
will be shown to add nothing to the
qualitative conclusions
based on 100 realizations.
There are several restrictions on the application of the correlation distortion method. The
static transformation
suggested in equation (4) is appropriate
for processes in which the
non-Gaussian
behavior can be adequately limited to a non-zero skewness and a kurtosis
not equal to three. For processes for which moments beyond fourth order are necessary
to adequately
describe the non-Gaussian
behavior, the method loses accuracy, as this
higher order information
is distorted through the inverse and forward transformations.
Further, the solution of equation (4) for R,,(z) is not guaranteed to be positive definite for all
L(r).
Direct transformation
An alternative
to the correlation-based
approach
Gaussian time history rather than its autocorrelation.
simulations of the sample process. The non-Gaussian
its Gaussian underlying form, u(x), through
U(X) =
[Jm
is to begin with a sample of a nonThe schematic in Fig. 3 then provides
sample process, x(t), is transformed to
+ i'(~)l"~ - [Jm
-
K$1"3
- a,
where
63
“=3fi,
c = (b - 1 - a2)3,
b=+
(6)
4
and the other parameters are defined after equation (4). Subsequently,
linear simulations
created through standard techniques based on the target spectrum of the Gaussian process,
G,,(W), are transformed
back to the non-Gaussian
parent form through eqn (5).
The shortcoming
of this direct transformation
technique is that the simulated
nonGaussian signal power spectrum does not match the sample non-Gaussian
spectrum to
a satisfactory
degree for the example sample processes we have used in subsequent
examples. This distortion may stem from the inability of the truncated Hermite moment
transformation
in eqn (6) to produce a Gaussian signal for cases when the parent signal is
highly non-Gaussian.
Specifically, non-Gaussian
behavior
requiring
moments
beyond
fourth order for characterization
are not addressed, and their presence distorts the static
transformation
from non-Gaussian
to Gaussian. The linear simulation
is then based on
a target spectrum derived from a process which is assumed Gaussian, but is not. It is at this
point, indicated in Fig. 3 by the dashed box, where the frequency information
is distorted,
and results in poor simulations.
One option for improving results is to add terms to the
Hermite series until a Gaussian transformation
is achieved. This may require a different
number of terms to achieve accuracy for varying input sample signals, and leads to a very
complex solution for U(X) in the higher order equivalent of equation (6).
An example of the potential for distortion using the direct transformation
method is seen
in Fig. 4. The top figures compare the same measured wind pressure signal seen in Fig.
2 with a direct transformation
simulation of the signal. The power spectral density and PDF
Class of non-normal
Fig. 3. Schematic
random
processes
of the direct transformation
605
method
Fig. 4. Measured wind pressure signal (top left), a direct transformation
simulation (top right), and
power spectral density and pdf of the measured data and ensemble of 100 simulations.
of the measured data are compared with an ensemble of 100 realizations in the lower figures
and show poor agreement. Clearly, this method is not applicable to this sample process.
Modi$ed direct transformation
A modification
is now suggested, shown in Fig. 5, to remove this distortion in the direct
transformation
method. Referring to equation
(6) it can be seen that the governing
parameters
L3, &, a, b, c, I, and thus u(x), are dependent
on the skewness and kurtosis,
y3 and y4. Since it is required that the process, u(x), be Gaussian in order to avoid the
distortion effects discussed above, y3 and y4 may be treated as adjustable input parameters
in order to force the transformed
process, u(x), to be Gaussian in terms of the third and
fourth moments. Optimization
of these two parameters is based on the minimization
of the
function
min(& + Y:.),
(7)
K. R. Gurley
606
Fig. 5. Schematic
-
of the modified
et al.
direct transformation
method
measured data zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFE
0
12
3
Fig. 6. Measured wind pressure signal (top left), a direct transformation
simulation (top right), and
power spectral density and pdf of the measured data and ensemble of 100 simulations.
where y3., y4_ are the skewness and kurtosis of the transformed process U(X). The optimized
input parameters y3 and y4 now provide a Gaussian process in terms of third and fourth
moments, and the linear simulations do not contain distortion of the frequency content. The
same parameters are used to transform back to a non-Gaussian
simulation whose PDF and
power spectral density closely match those of the sample process. This correction
is
essentially a quantification
of the error in truncating the Hermite series after the third term.
An example of the improvement
afforded by the modified direct transformation
method
is demonstrated
in Fig. 6. Again the measured pressure trace in the top left is simulated and
displayed in the top right. The power spectral density and PDF of the data and simulations
are shown in the bottom figures. Table 1 shows an improved abihty to match higher order
statistics compared with the correlation distortion method. By observing the positive and
Class of non-normal
I
0
I
500
random
I
I
1000
1500
processes
607
I
I zyxwvutsrqponmlkjih
3000
I
2000
2500
Modified Direct Transformation Simulation
I I
I
I
I
I
I
I
500
0
I
1000
1500
2000
2500
3000
2500
3000
Direct Transformation Simulation
0
500
Fig. 7. Measured
-
1000
TLP response,
modified
1500
2000
direct transformation
tions.
0
measured data
and direct transformation
measured data _________.
direct si&lation
modified direct simulation
0.6-
_________.direct simulation
modified direct simulation
14
0.5-
12
in
0.41
simula-
,’ ’
0.12
03
0
0.05
Fig. 8. Power
0.1
0.15
spectral
density
0.2
0.25
-3
-2
-1
and pdf of measured TLP response
realizations.
0
1
2
3
signal and ensemble
4
5
of 2000
negative extreme behavior, as well as fluctuation amplitude close to the mean, the modified
direct simulation can be seen to emulate the characteristics
of the measured process better
than correlation
distortion.
This behavior is quantified
by the kurtosis and standard
deviation, which match well with the data (Table 1).
A second example further demonstrates
the performance
of the modified direct transformation. In this case the sample process to be simulated is the measured response of
a model tension leg platform (TLP) under a random wind and wave field in a test facility.
The response is highly non-Gaussian
and has two dominant
frequencies. Figure 7 shows
a portion of the sample measurement,
a simulation
using modified direct transformation,
and a direct transformation
simulation in the top, middle, and bottom plots, respectively.
Figure 8 is a comparison
of the PDF and power spectral density of the sample and 2000
K. R. Gurley zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDC
et al
608
Table 2. Statistics
of measured
TLP response data and ensemble averaged
realizations, and (#) = 2000 realizations
Measured TLP data
Ensemble of 100 and (2000)
modified direct simulations
Ensemble of 100 and (2000)
direct transformation
simulations
simulated
data,
# = 100
Standard
deviation
Coefficient of
skewness
Coefficient
kurtosis
1.0
0.9720 (0.9690)
0.5165
0.8187 (0.8298)
3.7455
4.2127 (4.2650)
0.9633 (0.9672)
0.8419 (0.7546)
7.4672 (7.1469)
of
realizations
of the simulations.
Table 2 lists statistics from the data, an ensemble of 100
realizations,
and an ensemble of 2000 realizations
in parentheses. The direct transformation provides simulations
whose skewness characteristics
adequately
match the sample
(Table 2). However, large negative excursions in the realization
are not observed in the
sample, and lead to a significantly higher kurtosis (Table 2) as well as a poor fit of the PDF
to the data, most importantly
in the negative tail region. The modified direct transformation
provides realizations
which match the sample PDF well, particularly
in terms of positive
and negative extreme behavior.
The modified direct transformation
method, for the two sample processes considered, is
able to provide simulations which match the PDF and power spectral density of the sample
process, and match the scalar representations
of higher order statistics (skewness and
kurtosis) through fourth order. Later, an example will be presented where these comparisons are not as favorable.
The shortcoming
of any static transformation
is its inability to retain the phase interaction among related frequency
components.
The bispectrum
is a representation
of
the quadratically
phase coupled frequency components. Just as the power spectral density is
the distribution
of the variance of a signal with respect to frequency, the bispectrum is the
distribution
of skewness with respect to frequency pairs. Although the modified direct
transformation
is able to replicate the volume under the bispectrum, i.e. the skewness, it
is not able to correctly match the distribution
of skewness with respect to frequency.
Figure 9 compares the bispectrum contour of the sample TLP response process with that of
the direct and modified direct transformations.
Neither simulation
is able to adequately
match the shape of the sample bispectrum.
The non-linear
processes considered in this
study can be described by a quadratic
form [6,16,17].
Although higher-order
spectra
beyond the bispectrum may be calculated, it is only necessary to show the existence of the
bispectrum to demonstrate
the non-Gaussian
behavior of quadratic processes.
When the only available information
is a sample of the final process to be simulated (e.g.
wind pressure on a building face), this static transformation
method is quick and promising.
However, if more information
is available (e.g. the upwind wind velocity), it is possible to
better simulate the desired process by establishing a system identification
model between, in
this case, velocity and pressure. The limitation
of static transformation
techniques
is
overcome by the application
of memory-based
system identification
models.
TRANSFORMATIONS
WITH
MEMORY
When input and output are available, the non-linear
relationship
between current time
output and previous time input and output may be modeled through a variety of system
identification
techniques. A model that relates Gaussian input to non-Gaussian
output may
easily be applied to the simulation
of the output process by passing simulations
of the
Gaussian input through the model. The important
phase characteristics
of the output are
then retained in the simulation through the memory transformation,
assuming the model is
accurate. Here we consider the application
of a Volterra series formulation
for polynomial
non-linear
systems, as well as the empirical development
of discrete non-linear
differential
models through
NARMAX
(non-linear
autoregressive
moving average models with
exogenous inputs) and neural network models. Examples are provided using a two-term
Volterra series with analytical kernels, and an empirical neural network model.
Class of non-normal
random
processes
609
Measured non-Gaussian process
_
_
frequency
Modified direct transformation
frequency
Direct transformation
friquency
Fig. 9. Contours
of the bispectrum
of the measured TLP response and the bispectrum
of an
ensemble of 2000 simulations using the modified direct, and the direct transformation
methods.
System identification
models for non-Gaussian
input and non-Gaussian
output are
readily available, but are not easily adaptable for simulation
purposes. We wish to use
transformations
to relate easily attained Gaussian simulations to the desired non-Gaussian
process. The inability to easily simulate non-Gaussian
system input is addressed by the
hybrid application
of a static transformation
in combination
with a neural network. zyxwvutsrqponmlk
Volterra series model
In the Volterra series formulation,
the input-output
relationship
may be expressed in
terms of a hierarchy of linear, quadratic and higher-order
transfer functions or impulse
response functions [17-191. These transfer functions can be determined from experimental
data or from theoretical considerations.
For example, a non-linear
system modeled by
Volterra’s stochastic series expansion is described by
+
h,(r,,
z2,
~+(t
- z,)x(t
- ~~)x(t - z,)drI
dT2 drs + 3 . . ,
(8)
K. R. Gurley
610
where
sponse
The
in the
et al.
h,(r), hz(rl, r2) and h3(r1, r2, r3) are the first, second and third-order
impulse refunctions.
Fourier transform of the Volterra series expansion in equation (8) gives the response
frequency domain as
Y(h) = ffl(fi)X(h)
c
+
~2(flLf2)x(.fl)x(f2)
fl+f2=f;
+
c
~3(fl,f2,f3)X(fl)X(f2)xo
+ . ‘.
(9)
fl+f2+ls=f,
The Volterra series model in the frequency domain [equation
(9)] lends itself to the
simulation of non-linear
processes for which the Volterra kernels are available or may be
estimated. A non-Gaussian
signal resulting from a quadratic transformation
of a Gaussian
process may be simulated by the addition of second-order
contributions
to the complex
spectral amplitude components
at the appropriate
sum and difference frequencies before
inverse Fourier transforming
the sequence to the time domain. These second-order
contributions are formed from the products of pairs of linear Fourier components
with the
quadratic transfer function (QTF) in the frequency domain, and correlate the phase between
various frequency components
to a degree weighted by the QTF. The memory retained by
convolution
with the QTF facilitates the simulation of processes that are able to match not
only the power spectrum and PDF of the parent process, but the bispectrum as well, e.g.
Peinelt and Bucher [20].
The estimation of the higher-order
transfer functions in equation (9) requires the calculation of the cross-bispectrum
and cross-trispectrum
of the input and output processes. As
examples of the utility of bispectra and trispectra in frequency domain analyses, consider
two types of non-linear functions of a zero-mean, Gaussian random process, u(t). Functions,
g(u(t)), for which all odd-order moments vanish, will be considered statistically symmetric
non-linearities,
while those for which, in general, all moments are non-zero will be considered statistically asymmetric non-linearities.
For instance, g(u(t)) = n3(t) is a statistically
symmetric non-linearity,
whereas g(u(t)) = Us is statistically asymmetric. In the statistical
characterization
of the statistically asymmetric non-linearity,
we expect non-Gaussianity
in
the form of non-zero skewness, which we can characterize in terms of frequency pairs via the
bispectrum.
Indeed, for some cases, we can successfully employ a technique known as
equivalent statistical quadratization,
which retains memory by employing a Volterra series
approach in the frequency domain, to approximate
a more complicated statistically asymmetric non-linearity
as a quadratic
polynomial
for the determination
of higher-order
statistics [l&21]. On the other hand, for the case of a statistically symmetric non-linearity,
we expect the skewness to vanish and with it, the bispectrum. Hence, we must turn to the
trispectrum,
the tri-variate frequency domain representation
of the kurtosis, to gain any
higher-order
statistical information
about the non-linear
process. While the trispectrum
would supplement
the statistical evidence of the non-Gaussianity
of an asymmetric nonlinearity, it is not necessary in such a case as it is in the case of a symmetric non-linearity.
Again, under certain circumstances,
we may approximate
more complex symmetric nonlinearities in polynomial
forms containing
only linear and cubic terms via equivalent
statistical cubicization.
Implementing
the Volterra framework as described above, we can
approximate
higher-order
statistics in situations
involving symmetric non-linearities
as
well.
When the input x(n) and output y(n) of a system is available, the information
can be used
to estimate the Volterra kernels in equation (9) directly. The first and second order transfer
functions are given by
H
(f,)
1
=
I
<y(h)x*m>
(10)
(IXLfP>
and
1 (x*(fl)x*u2) Y(fl +f2))
ff2(flJ-2) = j
<IX(fl)X(f2)l”)
’
(11)
Class of non-normal
random
processes
611
where ( ) is the expected value operator. Equation (9) is then applied to simulate y(n) based
on linear simulations
of x(n). The formulation
for the QTF given in equation (11) is under
the assumption
of a Gaussian
input process X(f). The linear and quadratic
transfer
functions
can also be estimated
for a general random
input, i.e. without assuming
particular
statistics of the input [23325]. For systems where the governing
differential
equation is known, several methods of analytically
approximating
the Volterra kernels
are available, including variational
expansion, harmonic probing, and successive approximation [26].
Conceptually
the Volterra series may easily be extended to the simulation of non-linear
processes beyond second order [equation
(9)], although considerable
computation
time
is added by convolution
of the Fourier components
with transfer functions
beyond
quadratic.
Also, the acquisition
of higher-order
transfer functions from measured data
becomes difficult and the number of parameters
necessary to describe them becomes
prohibitive.
NARMAX
model
When the system is not yet described, system identification
techniques may be used to
approximate
the Volterra kernels directly, or to develop a discrete governing differential
equation. For the latter, NARMAX model algorithms
have been developed to identify
non-linear
systems [27]. The Volterra kernels may then be developed from the parametric
NARMAX model through harmonic probing, etc. [24].
The discrete differential equation provided by the NARMAX model may be used as the
stand-alone
representation
of the system. Convenient
algorithms
are available that use
input/output
data to define a polynomial
model. These algorithms identify the relevant
terms in an initial model consisting of all possible combinations
of input, output, and noise
terms up to a user-specified polynomial
order and maximum lag. For example, a model
specified as second-order
with two delays begins as
Y(n) =f(Y(n
- I), Y(n - 2), x(n - l), x(n - 2), e(n - l), e(n - 2)),
(12)
wheref( ) is the sum of all first and second order combinations
of the arguments.
Discussions
on NARMAX
theory are available in the literature, and a sampling of
practical and efficient algorithms may be found in Chen et al. [28], and Billings and Tsang
[24]. NARMAX
provides a more flexible representation
of a non-linear
system, and
generally requires fewer parameters than a Volterra model.
Neural network model
Another recently developed approach to non-linear
system identification
is the application of neural networks. A multi-layered
set of processing elements receives input information and uses the desired final output information
to adjust a weighting factor between each
of the elements.
Figure
10 shows such a network
with three weighting
layers
Wij(m), m=1...3,
where i=l...N,,
j=l...N,_r,
and N, and N,-r
are the
number of elements in the mth and the m - lth layers, respectively. The network in Fig. 10
has two hidden element layers ai
and a;(2) between the input and output layers ai
and
a;(3). In this example the input layer consists of the input occurring at the same time as the
current output from a,(3), and the two inputs preceding this lead input (a two delay input
system). Wi,Jm) then represents the weighting of the output from the element aj(m - 1)
before its input to element ai(
The output of each element is a non-linear function of the
weighted linear sum of the output from each of the elements in the previous layer as in
Kung [29]: zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
bi(nZ)= ‘2’
Wij(F?l)Llj(m- 1) + @i(m);
(13) zyxwvutsrqpo
1 di<N,;
(14)
j=l
ai(m)=f(bi(m))
l<m<3,
K. R. Gurley
612
et al
output layer
wij
C3)
hidden layer
wj,(2)
hidden layer
Wi,(l)
input layer
Fig. 10. Multilayer
neural network
0
15
Fig. 11. Measured
5
16
with three weighting layers and two hidden layers (adapted
Kung [29]).
15
10
17
TLP response
18
20
19
20
25
21
signal, trained and predicted
up of the training section.
22
40
35
30
23
neural network
from
24
output,
25
and a close
where 0,(m) is a threshold value fixed for each ai(
Various non-linear
functions may be
applied at the elements, under the restriction that the output must be limited to 0 df(bi)
d 1. One commonly applied function is the sigmoid function:
(15)
where (T is a parameter to control the shape Off(bi).
The element weights in the neural network are adjusted iteratively, commonly
with
a back propagation
scheme, which minimizes the error between the resulting and desired
final output. This is known as the training phase, in which the optimum model parameters zyxwvutsrqponm
Wij(?ll.), m = 1 . . . M are identified, where M is the number of network layers, and M = 3
for the example in Fig. 10 [29].
An example of the development
of a neural network is shown in Fig. 11. The measured
TLP response data in Fig. 7 is used as the input to a non-linear
difference equation, and
Class of non-normal
random
processes
613
a neural network is used to identify this input/output
system. The neural network was
trained on a 10 s span of input/output
from 15 to 25 s as identified in the top figure, and
shown alone in the bottom figure in Fig. 11. The actual desired system output and the
neural network estimate are both in the figures, and coincide almost exactly. The entire 40 s
input record is then passed through the model, and is shown in the top figure to predict the
actual output from the non-linear
difference equation extremely well. This accurate prediction capability will be used later for simulation purposes. zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQ
Applications
When the convenience
of having measured input and output is available, the Volterra,
NARMAX and neural network models discussed in the previous section may be used as
simulation tools when the input is Gaussian. The input is simulated and passed through the
prediction model to produce a non-Gaussian
simulation. This technique is much more time
consuming than simply applying a static transformation
technique to the output alone, but
has the advantage of memory built into the model.
A sample of a non-linear
simulation
using a Volterra series model is shown in Fig. 12.
This realization
is the surface elevation of gravity waves, with the non-Gaussian
train
showing the characteristic
high peaks and shallow troughs. In this case the second-order
Volterra kernel is analytically
derived [3&33] and referred to as a non-linear
interaction
matrix (NIM). The NIM relates a quadratic
non-Gaussian
process to its underlying
Gaussian
process. In terms of equation
(9) Y(f) is the desired non-Gaussian
wave
elevation, X(f) is the underlying linear sea state, and H,(f) is unity. X(f) is first simulated,
then used to generate the second-order contributions.
The matching of the realizations with
the desired target QTF is shown in Fig. 13, where the recovered QTF is an ensemble of 1000
realizations.
A non-linear
transformation
of Gaussian wave elevation is used for a neural network
example. The system input is a linear wave train simulated based on a JONSWAP spectrum
with a peakedness of 5, and a peak frequency of 0.05 Hz. The non-linear
output, F(n), is
generated from the linear wave train, q(n), by a generic non-linear
function:
F(n) = O.lq(n - 2) + 0.4$(n
- 2) + O.ly(n - 1) + 0.5y2(n - 1)
(16)
+ 0.2y(n) + 0.6$(n).
-2’
0
I
10
Fig. 12. Realization
I
20
of a Gaussian
I
30
I
40
I
50
I
60
I
70
and non-Gaussian
wave height generated
a non-linear interaction
matrix.
I
80
I
90
by Volterra
I
100
series using
614
K. R. Gurley
et u1
0.6,
0.4,
recovered QTF
Fig. 13. Comparison
Table
of target
3. Statistics
target QTF
QTF
applied in Fig. 12, and the recovered
realization ensemble.
of measured
non-linear
simulated
Measured wave data
Ensemble of 10 modified direct
transformation
simulations
Ensemble of 10 neural network
simulations
/
wave
data
process
and
QTF
ensemble
from a 1000
averaged
Standard
Coefficient of
skewness
0.4950
0.3948
2.2800
1.8911
9.8329
8.5284
0.4692
2.1256
8.6640
Sample Output
I
I
I
Coefficient
kurtosis
I
of
I
neural network Simulation
I
I
0
500
41
I
I
1000
I
1500
I
I
I
I
I
I
I
/
2000
2500
3000
3500
I
4000
I
I
Modified Direct Transformation Simulation
I
0
I
I
I
I
I
I
I
I
I
I
I
I
500
1000
1500
2000
2500
3000
3500
1
4000
Fig. 14. Sample output from Gaussian sea state input using equation 16 (top), a simulation using
a neural network trained on the sample input/output
(middle), and a simulation using modified
direct transformation.
Class of non-normal
bispectrum
of sample output
random
processes
bispectrum
615
of neural network
simulation
0.2
-0.21
-0.2
:
-0.1
0
0.1
I
0.2
freq. (Hz)
bispectrum
of modified direct simulation
0.21
-0.2
-0.2
I
:
-0.1
0
0.1
Fig. 15. Bispectrum contour of equation
realizations (top right), and bispectrum
bispectrum
0.2
16 output (top left), bispectrum contour of 10 neural network
contour of 10 modified direct transformation
realizations.
of sample output
bispectrum
of neural network
simulation
freq. (Hz)
bispectrum
of modified direct simulation
Fig. 16. Isometric view of Fig. 15. Bispectrum of equation 16 output (top left), bispcctrum of 10 neural
network realizations (top right), and bispectrum of 10 modified direct transformation
realizations.
A neural network with two delays is trained to model the input/output
from 4096 data
points. This model is then used to simulate realizations
of the output in equation (16) by
passing Gaussian simulations of the input, if, through it. The modified direct transformation
is also used to simulate the output directly, without knowledge of the input. A comparison
616
K. R. Gurley
et al
of statistical results is presented in Table 3, where it can be seen that the modified direct
transformation
does not match the statistics as well as in previous examples.
Figure 14 presents the original sample output in the top figure, and a neural network and
modified direct transformation
realization in the next two, respectively. The process being
considered is a quadratic transformation
of a Gaussian process. In order to demonstrate
its
non-Gaussian
nature it is sufficient to consider the bispectrum. Figure 15 shows a contour
representation
of the bispectrum of the sample output process, and of an ensemble average
of 10 realizations using the neural network and modified direct transformation
models. At
first glance, the modified direct transformation
simulation bispectrum contour appears only
slightly different from that of the neural network simulation
and the sample, which are
almost identical. However, the bispectrum
from the modified direct transformation
is
significantly different, as seen in an isometric view of the bispectra in Fig. 16. This figure also
shows the neural network bispectrum to closely resemble that of the original sample, due to
the memory retention.
CONCLUSIONS
A class of non-normal
processes are simulated based on information
from the sample of
a process. Static transformation
techniques
are applied to perform simulations
of the
underlying Gaussian time or autocorrelation
sample. An optimization
procedure is used to
overcome errors associated with the truncation
of static transformations.
Several examples
are presented to demonstrate
the utility of this method. The inability of static transforms to
retain the specific structure of the phase relationship
between frequency components
is
addressed by the application of memory models. A Volterra series up to second order with
analytical kernels is employed to simulate a non-Gaussian
sea state. A neural network
system identification
model is utilized for simulation
of output when system input is
Gaussian
wave elevation. This process is also simulated
by the modified direct static
transformation
method using only the output sample process. It is demonstrated
that the
memory model is better able to achieve the shape and magnitude of the bispectrum of the
original sample.
Acknowledyements-The
support for this work was provided in part by ONR grant N00014-93-1-0761,
and NSF
grants CMS9402196 and CMS95-03779. The first and third authors were uartiallv zyxwvutsrqponmlkjihgfedcbaZYXWVU
SuDDOrted
bv a Denartment
of
Education
GAANN Fellowship and a travel grant from the Institute of’Engineering'Mechanics, U‘niversity of
Innsbruck,
during this study. Neural network software was developed by Ioannis Konstantopoulos
under the
guidance of Drs Nicos Makris and P. J. Antsaklis. Their cooperation
is greatly appreciated.
REFERENCES
1. M. Shinozuka, Simulation of multivariate
and multidimensional
random processes. J. Acoust. Sot. Am. 49,
3577368 (1971).
2. M. P. Mignolet and P. D. Spanos, Recursive simulation of stationary multivariate
random processes. J. Appl. zyxwvutsrqpon
M ech. 54, 674- 687 (1987).
3. Y. Li and A. Kareem, Simulation of multivariate
random processes: hybrid DFT and digital filtering approach.
J. Enyny Mech. ASCE 119, 1078-1098 (1993).
4. T. T. Soong and M. Grigoriu,
Random Vibration of M echanical and Srructural Systems. Prentice-Hall,
Englewood Cliffs, NJ (1993).
5. M. K. Ochi, Non-Gaussian
random processes in ocean engineering. Prob. Engny M ech. 1, 28- 39 (1986).
6. A. Kareem, K. Gurley and M. Tognarelli,
Advanced analysis and simulation
tools for wind engineering.
International Association for W ind Eny ineering, Proc. Ninth Inc. Conf: W ind Engineering, Vol. 5. Wiley
Eastern, New Delhi (1995).
7. M. Grigoriu, Crossing of non-Gaussian
translation
process. J. Engng Mech. ASCE 110(4), 610-620 (1984).
8. S. R. Winterstein,
Nonlinear
vibration models for extremes and fatigue. J. Engng Mech. ASCE 114(10),
1772-1790 (1988).
9. R. N. Iyengar and 0. R. Jaiswal, A new model for non-Gaussian
random excitations. Prob. Enyny Mech. 8,
281-287 (1993).
10. R. Deutsch, Nonlinear Transformations of Random Processes. Prentice-Hall,
Englewood Cliffs, NJ (1962).
11. D. Ammon, Approximation
and generation of Gaussian and non-Gaussian
stationary processes. Struct. Safety
8, 153-160 (1990).
12. F. Yamazaki and M. Shinozuka, Digital generation of non-Gaussian
stochastic fields. J. Engny Mech. ASCE
114(7), 118331197 (1988).
13. M. Grigoriu, Applied non- Gaussian Processes. Prentice-Hall,
Englewood Cliffs, NJ (1995).
14. D. A. Conner and J. L. Hammond,
Modelling of stochastic system inputs having prescribed distribution
and
covariance functions. Appl. M ath. M odel/. 3(2), 67- 69 (1979).
15. G. E. Johnson, Constructions
of particular random process. Proc. IEEE 82(2), 27&285 (1994).
Class of non-normal
16. G. I. Schueller
Structural
and C. G. Bucher, Non-Gaussian
Dynamics,
Progress
in Theory
random
response
processes
617
of systems under dynamic excitation. In zyxwvutsrqponmlkjihg
Stochastic
Ariaratnam,
Schueller and Elishakoff (eds),
and Applications,
pp. 219-239. Elsevier Applied Science, London (1988).
17. A. Kareem and Y. Li, On modelling the nonlinear relationship
between random fields by means of higherorder spectra. In Probabilistic M ethods in Civil Engineering, P. D. Spanos (ed.), pp. 384-387. ASCE, New York
(1988).
18. P. D. Spanos and M. G. Donley, Equivalent statistical quadratization
for nonlinear systems. J. Engng Mech.
ASCE 117(6), 1289-1309 (1991).
19. M. Schetzen, The Volterra and zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA
Wiener Theories of NonlinearSystems. Wiley, New York (1980).
20. R. H. Peinelt and C. G. Bucher, Spectral analysis and synthesis of non-Gaussian
processes. In Structural Safety
and Reliability, G. I. Schueller, M. Shinozuka
and J. P. Yao (eds), pp. 195-200. Balkema, Rotterdam
(1994).
21. J. Zhao and A. Kareem, Response statistics of tension leg platforms under wind and wave loads: a statistical
quadratization
approach. ICOSSAR, Austria (1993).
22. A. Kareem and J. Zhao, Stochastic response analysis of tension leg platforms: a statistical auadratization
and
cubicization
approach. In Proc OMAN ‘94 Conj!, Vol. 1. ASMET New York (1994).
_
23. S. W. Nam, E. J. Powers and S. B. Kim, Applications
of digital polyspectral
analysis of non-linear system
identification.
In Proc. 2nd IASTED Int. Svmp. Sicmal Processinu and its Aonlications. Gold Coast. Australia.
pp. 133-136 (1990).
24. S. A. Billings and K. M. Tsang, Spectral analysis for non-linear systems, Part I: parametric non-linear spectral
analysis. Mech. Systems Signal Process. 3(4), 319-339 (1989).
25. J. S.Bendat and A. G. Piersol, Random Data Analysis dnd M easurement Procedures. Wiley, New York (1986).
26. M. Wright and J. K. Hammond,
The convergence
of volterra series solutions to nonlinear
differential
equations. In Structural Dynamics: Recent Advances. Proc. 4th Inc. Conf., M. Petyt, H. F. Wolfe and C. Mei
(eds), pp. 422-43 1. Elsevier Applied Science, London (1990).
27. I. J. Leontaritis and S. A. Billings, Input-Output
parametric models for non-linear systems, Part I: deterministic non-linear systems. Inc. J. Control 41(2), 303-328 (1985).
28. S. Chen, S. A. Billings and W. Luo, Orthogonal
least squares methods and their application
to non-linear
system identification.
Int. J. Control 50(5), 1873-1896 (1989).
29. S. Y. Kung, Digital Neural Networks. Prentice-Hall,
Englewood Cliffs, NJ (1993).
30. R. T. Hudspeth and M. C. Chen, Digital simulation of nonlinear random waves. J. W aterways Port Coastal
Ocean Div. ASCE 105, 67-85 (1979).
31. K. Hasselmann,
On the nonlinear energy transfer in a gravity wave spectrum, Part I. /. Fluid Mech. 12,
481-500 (1962).
32. L. J. Tick, Nonlinear probability
models of ocean waves. In Ocean W ave Spectra, pp. 163-169. Prentice-Hall,
Englewood Cliffs, NJ (1963).
33. A. Kareem and Hsieh, Probabilistic
dynamic response of offshore platforms to wave loads. Technical Report
no. NDCEYI-I, Department
of Civil Engineering,
University of Notre Dame (1991).