Representation for the Slepian process with applications

Journal of Sound and Vibration (1988) 124(3),435-441

REPRESENTATION FOR THE SLEPIAN PROCESS WITH APPLICATIONS A. M. HAsoFER Department of Statistics, University of New South Wales, Kensington, New South Wales, Australia (Received 2J February 1987, and in revised form 12 November 1987)

An alternative representation for the Slepian process based on the principle of conditional simulation is obtained. It is particularly suitable for efficient simulation and is applied to the evaluation of excursion times and reliability functions. 1. INTRODUCTION

The Slepian model process is named after D. Slepian who introduced it [1]. It is a random function representation of the behaviour of a Gaussian process conditioned on events defined by its level crossings. Its explicit structure makes it well suited for probabilistic manipulations, finite approximations and asymptotic expansion. It is also used to derive optimal level crossing predictors and to study wave-characteristic distributions in random waves, with applications to fatique [2]. As explained at length in Slepian's original paper, it is the natural tool to use in studying the behaviour of a stationary Gaussian process after a level crossing, whether at the same level or at a different one. However, even though the Slepian process has an explicit form, it does not yield exact formulae for waiting time distributions, only approximations and bounds. To obtain exact results, resort must be had to simulation. Unfortunately, as usually expressed, the Slepian process involves a non-stationary component that is difficult to simulate precisely. In this paper, an alternative form of the Slepian process is given that involves just three random variables in addition to the original process. Moreover, the new representation permits an easy evaluation of the error induced in the Slepian process when only an approximate simulation of the original process is available. The new representation is applied to the evaluation of the excursion time and transient reliability function of the linear oscillator driven by white noise. It might be argued that ordinary simulation of the stationary underlying process is sufficient to solve the problems listed above since it is sufficient to simulate until, say, a level crossing occurs, and then study the behaviour of the process. There are two objections to this: (a) The simulation effort until a level crossing occurs is wasted; this is particularly acute if a high-level crossing is being studied; (b) sometimes the conditioning event has zero probability of occurring in any finite interval; This is, for example, the case when a linear oscillator is started from rest (see section 8); in that case the conditioning event g( t) = g'( t) =: 0 is unlikely ever to occur, no matter how long the process is simulated. The method advocated in this paper readily generalizes to conditioning events more complex than just level crossings, with the structure of the algorithm unchanged. 2. CONDITIONAL SIMULATION

Let X be a Gaussian vector with mean J.l..x and covariance matrix .l'11> and let Y be another Gaussian vector with mean J.l.. y and covariance matrix 1:22 , Let the cross-covariance 435 0022-460X/88/150435+ 07 $03.0010

© 1988 Academic Press Limited

436

A. M. HASOFER

matrix between Y and X be .r 21. Then it is well known (e.g., see reference [3]) that the conditional mean and variance of Y given that X = Xo are given by (1)

and (2)

where II2 = III' Let, for brevity, A = I21.r~IJ. It is required to simulate a Gaussian vector Y c that has the same distribution as Y given X = Xo. Proposition: To generate Y, proceed as follows: (1) generate vectors Y*, X* from the joint distribution of X, Y; (2) set (3)

Yr=Y*+A(Xo-X*).

Proof of the proposition. Clearly, Yc is Gaussian since it is a linear combination of two jointly Gaussian vectors. It remains to verify that Y c has the required mean and covariance matrix. Indeed

Cov (Yol = I

J.A.r = E(Y c) = J.A.y + A(Xo-J.A.x),

(4)

T

(5)

22 -

AI 12 - I

2 1A

+ AII1A T = .l'22-.l'21.r~11 I12,

as required. The principle of conditional simulation is to be found in special form in the paper by Journel [4].

3. THE SLEPIAN MODEL PROCESS

Let X( t), t E R, be a real-valued stationary Gaussian process with zero mean and covariance function ret), assumed twice differentiable. The Slepian process X,,(t) associated with X(t) is defined by ur(t) Zr'(t) ( ) X () t =-----+K t "

>"0

A2

(6) '

where Ao = reO) and >"2 = -r"(O) are the spectral moments of X (t) of order zero and two respectively, Z is a random variable having the Rayleigh density function Iz(z)=(z/>"2)e- z 'f2>.."

z;;;':O,

(7)

and K(t) is a zero-mean, non-stationary Gaussian stochastic process independent of Z, with covariance function

C ( t, u ) := r ( t - u )

r(t)r(u)

r'(t)r'(u)

An

A2

.

(8)

Roughly speaking, the process X,,(t) has the same probability law as {X(t)IX(O)=u,

X'(O)

~O}.

For a complete discussion of the exact meaning of the conditioning of X (t) as well as the exact conditions under which the above results hold see the book by Leadbetter et al. [5J. Suffice it to mention that equation (6) is obtained by noticing that

ur( r) zr' (t) E{X(t)IX(O)= u,X'(O)= z} : = - - >"0

A2

(9)

437

SLEPIAN PROCESS WITH APPLICATIONS

and that COy

{X(t), X(U) \X(O)

= u, X'(O) = Z} = C(/, u)

(10)

independently of z, and then deconditioning with respect to X'(O). Note now that simulation of Xu(t) requires the simulation of K(/), a new non-stationary process with a rather awkward convariance function C(/, u). In the next section it is shown that Xu(t) can actually be simulated by simply adding to the original process X(t) just three random variables, an easy and efficient algorithm. 4. SIMULATION OF THE SLEPIAN MODEL PROCESS BY CONDITIONAL SIMULATION OF X(t)

Suppose that there is available some method for simulating the original process X( I) by X"'(t), and that it is required to simulate the associated Slepian process XlI(t). The first step is to simulate XlI(t) conditional on X~(O) = z. For this purpose the proposition of section 2 is used. Set X = {X"'(O), X*'(O)}, X o = {u, z} and Y* = {X*( r), I> OJ. Upon noticing that -1

4- 1 11 -

(

Ao 0

and that E{X*(/)X*(O)} = r(t), E{X*(t)X*'(O)} = -r'(t) the following representation is obtained:

xtz(t) = X*( t) + [r( t)/ '\0][ u - X*(O)] - [r'( t)/ A2 ] [ z- X*'(O)].

(11)

Here xtz(t) is a simulation of X~(t) conditional on X~'(O) = z. This can be rewritten in the form

X~z(t) = U~~) _ zr~~t) + [ X*(t) -

X*(O)

r~:) + X*'(O) r';2 J.

(12)

which is exactly of the same form as equation (6) with K(/) replaced by the expression in square brackets. It is in fact easily verified that the latter expression has zero mean and covariance function C(t, u) as in equation (8). Upon deconditioning with respect to z, the following expression is finally obtained for X~(t):

X~(t) =

ur(t) _ Zr'(t) + [X*(t) _ X*(O) r(t) + X*'(O) r'(t)J. Ao

A2

AO

(13)

,\ 2

Here, as before, Z is a Rayleigh distributed random variable, independent of X (t) and having the density function given in equation (7). 5. APPROXIMATE SIMULATION OF THE SLEPIAN MODEL PROCESS

Denote now the process of interest by ( t), with a simulation (*( t), and suppose that (*(/) can be represented in the form e(/) = X*(/)+ 1')*(t),

(14)

where X*(t) is a zero-mean Gaussian stationary process. It is assumed that X*(/) can be easily simulated, and that 1') *( t) is a residual process that cannot be easily simulated. However, bounds on the size of 1')*(1) can be calculated. For situations where such a problem arises see the papers by Hasofer and Ghahreman [6] and by Hasofer [7]. Let

438

A. M. HASOFER

r( t) be the covariance function of g( t), assumed known, and Slepian process associated with ~(t), namely

~~( r)

a simulation of the

~~= ur(t) _ Zr'(t) +[e(t)-g*(O) ret) +r'(O) r'(t)]. Ao

A2

Ao

(15)

A2

Set

x~*( t) = ur( t) _ Zr'( t) + [x*( t) _ X*(O) Ao

Az

r( t) + X*'(O) r'( t)]. Ao A2

(16)

It is to be noted that X~*(O) = u, and X~*'(O) = Z. Then

ret)

r'(t)

~~( t) - X~*(t) = 11*( t) -11*(0) ~+ 7]*'(0) ~ .

(17)

Suppose that the time interval of interest is (0, T). Then max 1~~(t)-x~*(t)lo=; max 17]*(t)I+I7]*(O)I+I7]*'(O)I(A o/ Az) l/ z, O'Ei,""T

O"',""T

(18)

since Ir(t)1 ~ Ao, and Ir'(t)1 0=; (AoA 2 ) 1/ 2, by Schwarz's inequality. Let (T~, lT~ be the spectral moments of 'l']*(t) of order zero and two respectively. Then it is a standard result that P( max 0"'0;; T

1'l']*(t)I"""x)0=;2{1-.p(~)}+ T(T2e-x2/2(T~, lTo

(19)

71'CTo

where the first term is the instantaneous probability and the second term is the mean number of outcrossings in time T Finally the probability that the right side of equation (18) is larger than a given value can be bounded by repeatedly using the inequality P(X""" x+ k) - p(1 yl;;,: k)"';; P(X + y;;,: x)",;; P(X;;': x-k) + p(1 Y!;;,: k),

(20)

which holds for any random variables X, Y and any constants k;;,: 0, x [8]. 6. APPLICATION TO THE CALCULATION OF EXCURSION TIME

Let Wu be the excursion time at level u of the zero mean, stationary Gaussian process that is, the time between an upcrossing of level u and the next downcrossing. Let ~~(t) be the simulated Slepian process associated with get), and let ~(t):

Mg(t) = min g~(t). O~T~'

(21)

Then clearly (22)

Let now MAt) = min xt*(t)

(23)

OO!Ei'TEr

and (24)

Then clearly M,,(t) - MR(t)",;; Mg(t)

0=; Mx(t) +MR(t),

(25)

so that by using equations (18) and (20) upper and lower bounds for P( l¥.,;a.: t) can be calculated in terms of the excursion time probabilities of X**( t).

439

SLEPIAN PROCESS WITH APPLICATIONS

7. EXCURSION TIME OF THE LINEAR OSCILLATOR RESPONSE

The above method was applied to the stationary response of a linear oscillator driven by white noise, described by the equation f'(t)+2{wof(t)+w~g(t)= W(t).

(26)

It is well known that if the spectral density of the driving Gaussian white noise is taken

to be

JAw) = 1/21T, -00 < w < +ro,

(27)

then the spectral density of the stationary response X(t) is given by [9] J.(w) = (l/21T ){1/[ (w~ - W2)2+ 4{2W~W2]}.

(28)

The problem of the distribution of the level crossings of Ht), for various initial conditions, has been widely studied, and has an extensive literature. A glimpse of the current interest in the problem can be obtained from two recent papers: those of Spanos et al. [10] and Krenk et al. [11]. If the process is conditioned on an upcrossing at u, then the time to the first downcrossing is the excursion time. However, despite the large amount of effort invested, no explicit formula for the excursion time distribution been found, only bounds, approximations and asymptotic results, so that for exact results it is necessary to rely on simulation. Now it is true that the linear oscillator response can be simulated exactly at equispaced intervals by mean of an autoregressive moving average process (see the paper by Gersch and Liu [12]). But that approach to simulation does not yield any information on the fluctuation of the response between the equispaced points. Moreover, application of the simulation algorithm introduced above requires the calculation of the derivative of the response process at the origin. The value of that derivative cannot be obtained with known precision from the ARMA process. Simulation of stationary Gaussian processes by means of finite trigonometric polynomials has been advocated by Hasofer [7] and applied to the linear oscillator response by Hasofer and Ghahreman [6]. It has two advantages over the ARM A model: (a) it provides a continuous-time formula for the simulation, thus enabling the derivative to be calculated; (b) it provides an upper bound for the residuals of the process and its derivative. It has been shown by Hasofer and Ghahreman [6] that g(t) can be closely approximated by a random trigonometric polynomial of the form X(t)=

N, L n=N,

( nsrt n1Tt) an Vncos-+Wnsin2T 2T

(29)

over the interval (- T, T), where {Vn } and { Wn } are independent standard normal variables. For the values T=2, wo= 151T, and {= 1/15, the choice N, = 0, N 2 = 511 was made. The largest coefficients are in the range n = 54 to 64 and the largest is 2·85 x 10-6 • The original process g(t) can be written as (30)

g(t) = X(t)+ 1/(t),

with 'Y/(t) independent of X(t). For 'Y/(t) it was found that (J'~=2'OX 10- and (J'~= 1·80 x 10- 6 as against Ao = 1·126 X 10- 4 and A2 = 0-080. The distribution of the excursion time was estimated from 3000 simulations of X(t). Some of the values obtained together with their estimated standard deviations are given in Table 1. 9

440

A. M. HASOFER TABLE

1

Distribution of the excursion time W;, of X (t) at u =0,001 (ujJA o = O' 106 )

0·15 0·30 0·45 0·60 0·75 0·90 1'05 1·20 1·35 1·50 1·65 1-80

P(W
Standard deviation

0·008 0·022 0·055 0'109 0·241 0·479 0·803 0·924 0'964 0·981 0·989 0·994

0·0016 0·0027 0·0042 0·0057 0·0078 0·0091 0·0073 0·0048 0·0034 0·0025 0·0018 0'0014

Correction for the residual process turned out to be negligible. It is of interest to note that in this case the Slepian representation (11) turns out to have a very simple interpretation. Let w~ =: w~( 1 - 52), and solve equation (26) as a deterministic equation. The general solution ga(t) is given by to(t) = e -{Wo'[K) cos Wdt + K 2 sin Wdt] + 5s(t),

(31)

where ts(t) is the steady state solution and K 1 , K 2 are arbitrary constants. Upon introducing the initial conditions to(O) = u and to(O) = z, and using the expression R(t) = (lj4wdW~) e-'Wo'{(WdUWO) cos Wd!}

(I;;;'

0)

(32)

for the stationary covariance function, the solution becomes (33)

Thus the last two terms of equation (11) turn out to be simply the "transient" terms generated by the initial conditions and the excitation generating X*(t) before 1 = O. This interpretation is quite general since those two terms actually tend to zero as t tends to infinity, if r(t), r'(t) tend to zero as t tends to infinity. In the present case R(t) = r( t) since here the vector process {get), t'(t)} is Markovian, due to the white noise nature of the driving function Wet). 8. TRANSIENT RELIABILITY OF THE LINEAR OSCILLATOR

Suppose the linear oscillator described by equation (26) is started from rest: i.e., with g(O) = g'(O) = O. The reliability function L g ( I, u) is defined as the probability that {( t) remains within the band (-u, u) in (0, t). Evaluation of L~Ct, u) has been attempted by Corotis, Vanmarcke and Cornell [13], and more recently by Petocz [14] for the oscillator with ~ = 0'1 and u = (T(30/ wo), where (T(t) is the standard deviation of g( t) at t. Comparable results were calculated, by using the method of section 7, modified for the initial values tCO) = g'(O) = 0, together with the standard deviation of L~, from 3000 simulations. Results are given in Table 2. The results show that previous results tend to be too conservative. A step-by-step description of the calculations in sections 7 and 8 is available in reference [15].

SLEPIAN PROCESS WITH APPLICATIONS

441

TABLE 2

Comparison of reliability function L s (t, u) obtained by various methods t

CVC

P

H

S.D.

15 20 25 30

0·72 0·58 0·47 0·40

0·71 0·57 0·48 0·37

0·76 0·69 0'56 0·50

0·008 0·008 0·009 0·009

Wa

eve = Corotis et al.; P = Petocz; H = present paper; S.D. = standard deviation of present paper results. REFERENCES

1. D. SLEPIAN 1963 in Time Series Analysis (Editor M. Rosenblatt), 104-115. New York: John Wiley. On the zeros of Gaussian noise. 2. G. LINDGREN 1984 in Statistical Extremes and Applications (Editor J. Tiago de Oliveira), 261-285. Dordrecht, Holland: D. Reidel. Use and structure of Slepian model process for prediction and detection in crossing and extreme value theory. 3. C. R. RAO 1973 Linear Statistical Inference and its Applications 2. New York: John Wiley. 4. A. G. JOURNEL 1974 Economic Geology 69, 673-682. Geostatistics for conditional simulation of ore bodies. 5. M. R. LEADBETTER, G. LINDGREN and H. ROOTZEN 1983 Extremes and Related Properties of Random Sequences and Processes. New York: Springer-Verlag. 6. A. M. HASOFER and S. GHAHREMAN 1985 Proceedings of International Conference on Education, Practice and Promotion of Computational Methods in Engineering using Small Computers (EPMESC), Macao, 449-457. Simulation of random vibrations on a small computer. 7. A. M. HASOFER 1987 Journal ofSound and Vibration 112,283-293. Distribution of the maximum of a Gaussian process by a Monte Carlo method. 8. A. M. HASOFER 1982 in Essays in Statistical Science, Journal of Applied Probability Special Volume 19A, 333-344. Simple trigonometric models for narrow-band stationary processes. 9. Y. K. LIN 1967 Probabilistic Theory of Structural Dynamics. New York: McGraw-Hill. 10. P. T. D. SPANOS and G. P. SOLOMOS 1984 Journal of Engineering Mechanics, American Society Civil Engineers 110, 20-36. Barrier crossing due to transient excitation. 11. S. KRENK, H. D. MADSEN and P. H. MADSEN 1983 Journal of Engineering Mechanics, American Society of Civil Engineers, 109, 263-278. Stationary and transient envelope statistics. 12. W. GERSCH and R. S. Z. LIU 1976 Journal of Applied Mechanics 43, 159-165. Time series methods for the synthesis of random vibration systems. 13. R. B. COROTlS, E. H. VANMARCKE and C. A. CORNELL 1972 Journal af Bngineering Mechanics, American Society of Civil Engineers 98, 401-414. First passage of non-stationary random processes. 14. P. PETOCZ, 1981 Ph.D. Dissertation, University of New South Wales, Sydney, Australia. Upcrossings by oscillatory processes and their envelopes. 15. S. GHAHREMAN 1986 Ph.D. Dissertation, University of New South Wales, Sydney, Australia. Extrema of Gaussian process by simulation.