Second order stochastic simulation with specified correlation

Advances in Engineering Software 30 (1999) 489–494

Second order stochastic simulation with specified correlation D.B. Parkinson* Department of Engineering, Mechanical Engineering, University of Liverpool, Brownlow Hill, Liverpool, L69 3GH, UK Received 1 April 1996; received in revised form 1 October 1998; accepted 1 December 1998

Abstract A previously reported procedure for the simulation of normal stochastic processes using a combination of linear segments is extended by the use of non-linear second order curves, for which the second derivative process has known correlation function. The correlation function of the simulated process is deduced as a function of specified parameters, which are then determined to fit (by least squares optimisation) a given correlation function. Accuracy of fit is improved in comparison with the linear simulation. The sum of the resulting quadratic curves through sets of three random points is then the simulated process. Individual points in the process can be simulated without the need to store the whole process and as process parameters are determined prior to simulation the procedure is fast. These features render it suitable for use on microcomputers and for the simulation of large numbers of processes, for example, to simulate outcrossing statistics for vector processes in reliability and quality and process control theory. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Stochastic process; Quadratic process; Simulation

1. Introduction The simulation of random processes is used in a number of disciplines, for example, ocean wave simulation [1], simulation of earthquake forces or other naturally occurring phenomena, relevant to the safety of man made structures [2,3,5] to describe geological data [4] hydraulics and hydrology [6] and in the simulation of discrete event systems [7]. Simulations of this type are particularly useful in the assessment of theoretical predictions of rare, timedependent events (though the independent variables may also represent spatial location [4]). In the present case, the original motivating interest was the study of approximate failure risks in multivariate uncertainty models (reliability models), [8–11]. The latter references will give the reader numerous examples of the use of random processes in this context. Such models can, in principle, be extended to other applications [12–14]. Stochastic simulation has also found recent applications in the fields of biology and econometrics [15,16]. A further, possible application is in statistical signal processing. As most types of Stochastic process may be derived by transformation [17,18] from normal (Gaussian) processes and some non-stationary processes simulated by modification of an equivalent stationary process [2], the present

* Tel.: ⫹ 44-0151-794-4855; fax: ⫹ 44-0151-794-4848.

study is concerned with the improved determination of stationary normal processes. In reliability models, the requirement is usually to determine the probability of a vector Stochastic process passing out of a given safe region of the design space in a given time period [8,10]. Theoretical estimates of such a probability may be tested by simulation, constructing the vector process from a set of single variable processes [19]. As the event probabilities of interest are generally small and a large sample of outcrossing events is required for a statistical estimate of the required probability, long simulation runs are necessary. A number of existing methods [2,4,20,21] suffer the disadvantage in this context, of the need to store the simulated points together with long run times, where the simulation process is complex (see also Ref. [26]). To overcome these problems, a simple broken line process was derived in Ref. [19] and some improvements suggested in Ref. [22], where the possibility of extending the method to use higher order curves to replace the simple line segments was discussed. The original method produced a simulated process as a sum of a set of random processes consisting of linear segments with random end point values and incorporating a random shift for each component process [19,22]. The present article aims to improve on this method by replacing linear segments with second order curves with the aim of modelling a larger range of processes as defined by their (given) correlation functions to produce an improved

0965-9978/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0965-997 8(99)00003-4

490

D.B. Parkinson / Advances in Engineering Software 30 (1999) 489–494

simulated process curve. Using this method, the parameters of the simulation are determined prior to simulation and the latter then consists only in the generation of the process segments and their addition to form the final process. A fast simulation is therefore obtained and storage required is confined to temporary storage of the current segments parameters (see Section 4). There is no requirement to store the whole of the simulated process, which is a useful property when modelling rare event phenomena requiring many long simulation runs. It is also possible, by this method, to simulate selected regions or points of the process.

0–a :

Rb 00 …t† ˆ 1;

a–2a :

Rb 00 …t† ˆ

11 5 t ⫺ ; 6 6 a

2a–3a :

Rb 00 …t† ˆ

1 ; 6

3a–4a :

Rb 00 …t† ˆ

2 1 t ⫺ ; 3 6 a

…7†

To begin the derivation of the correlation function of the original process b it is noted that [25]:

2. Theory

E…b 00 …0†2 †Rb 00 …t† ˆ ⫺E…b 0 …0†2 †

Consider the stochastic process: 1 …wj⫹1 ⫺ 2wj ⫹ wj⫺1 †‰t ⫺ …j ⫺ 1†aŠ2 2a2

b…t ⫺ ka† ˆ

⫹

1 …⫺wj⫹1 ⫹ 4wj ⫺ 3wj⫺1 †‰t ⫺ …j ⫺ 1†aŠ 2a

E…b 0 …0†2 † ˆ ⫺

…1† obtained by fitting a quadratic through the three random points wj⫺1, wj, wj⫹1 separated in time (t) by interval a. k is a random shift parameter having values in the range zero to one and uniformally distributed. The wj are normal random variables having zero mean values and variances s 2, mutually independent and independent of k. The first and second derivative processes are:

b 00 …t ⫺ ka† ˆ

1 …⫺wj⫹1 ⫹ 4wj ⫺ 3wj⫺1 †: a

1 …wj⫹1 ⫺ 2wj ⫹ wj⫺1 †; a2

…2† …3†

For the latter process the covariance is:        6 s2 t s2 t t ⫹ 4 Fk 4 ⫺ ⫺ Fk 2 ⫺ ; Cb 00 ˆ 4 Fk 2 ⫺ a a a a a …4† where, Fk is the cumulative probability distribution of k. The variance is: 6 s2 ; a4

and the correlation function is then:     5 t 1 t ⫹ Fk 4 ⫺ : Rb 00 …t† ˆ Fk 2 ⫺ 6 a 6 a

6 s2 a4 R 00b 0 …0†

…9†

d 2 Rb 0 ˆ R 00b 0 …0†Rb 00 …t†; dt2

…10†

with Rb 00 …t† given in the aforementioned Eq. (7). Integration of (10), using Rb 0 …0† ˆ 1; …dRb 0 =dt† ˆ R 0b 0 …0† at t ˆ 0 and continuity of …dRb 0 =dt† and Rb 0 at a, 2a and 3a gives, for the specified ranges of t: Rb 0 …t† ˆ R 00b 0 …0†

0–a :

t2 ⫹ R 0b 0 …0†t ⫹ 1; 2

a–2a :

1 …wj⫹1 ⫺ 2wj ⫹ wj⫺1 †‰t ⫺ …j ⫺ 1†aŠ a2 ⫹

…8†

the variance of b 0 , and Rb 0 is given by:

…j ⫺ 1†a ⱕ t ⱕ …j ⫹ 1†a

b 0 …t ⫺ ka† ˆ

d2 Rb 0 …t† dt2

and as Rb 00 …0† ˆ 1 :

⫹ wj⫺1 for

E‰b 00 …0†2 Š ˆ

As Fk (z) ˆ z for 0 ⱕ z ⱕ 1, Fk (z) ˆ 0 for z ⬍ 0 and Fk (z) ˆ 1, for z ⱖ 1, Eq. (6) gives, for the specified ranges of t:

…5†

! 1 00 5a2 5 11 2 5 3 0 Rb 0 …t† ˆ R b 0 …0† ⫺ at ⫹ t ⫺ t R b 0 …0†t ⫹ 1; 6 6 2 2 6a

2a–3a : 1 35 15 t2 at ⫹ Rb 0 …t† ˆ R 00b 0 …0† ⫺ a2 ⫹ 2 6 6 2

! ⫹ R 0b 0 …0†t ⫹ 1;

3a–4a : 1 4 t3 Rb 0 …t† ˆ R 00b 0 …0† ⫺ a2 ⫹ 3at ⫹ 2t2 ⫺ 6 3 6a

! ⫹ R 0b 0 …0†t

⫹ 1: …11† Setting Rb 0 …4a† ˆ 0; then gives:

…6†

R 00b 0 …0† ˆ

⫺3 …1 ⫹ 4aR 0b 0 …0††: 16a2

…12†

D.B. Parkinson / Advances in Engineering Software 30 (1999) 489–494

491

Fig. 1. Correlation functions and process sample for example 1.

For t between a and 2a:

To obtain the required correlation function Rb(t), as d2 R E…b …0† †Rb 0 …t† ˆ ⫺E…b …0† † 2b ; dt 0

0

2

…13†

6 s2 ; a4 R 00b 0 …0†R 00b 0 …0†

…14†

2

and Rb 0 …0† ˆ 1 : E…b…0†2 † ˆ

For t between 2a and 3a: " # K 31 4 25 3 35 2 2 5 3 t4 ⫺ a ⫹ a t⫺ a t ⫹ at ⫹ Rb …t† ˆ 24 6 24 8 12 4

and d2 R b ˆ R 00b …0†Rb 0 …t† dt2

…15†

giving the variance of b(t) and a differential equation for Rb(t) which may be integrated using Rb(0) ˆ 1 and continuity of Rb(t) at a, 2a and 3a, yielding the following: For t between 0 and a: Rb …t† ˆ

Kt4 Jt3 t2 ⫹ ⫹ R 00b …0† ⫹ R 0b …0†t ⫹ 1; 23 6 2

where K ˆ R 00b …0†R 00b 0 …0† and J ˆ R 00b …0†R 0b 0 …0†:

"

a4 5a3 t 5a2 t2 5at3 11t4 ⫺ ⫺ ⫹ ⫹ 24 12 12 24 24 # 5t5 Jt3 t2 ⫹ ⫹ R 00b …0† ⫹ R 0b …0†t ⫹ 1: ⫺ 120a 6 2

K Rb …t† ˆ 6

⫹

Jt3 t2 ⫹ R 00b …0† ⫹ R 0b …0†t ⫹ 1; 6 2

and for t between 3a and 4a: K Rb …t† ˆ 6

"

11 4 a ⫺ 15 # t5 ⫹ ⫺ 120a

a3 t 2 1 1 ⫺ a2 t2 ⫹ at3 ⫹ t4 4 3 2 6 Jt3 t2 ⫹ R 00b …0† ⫹ R 0b …0†t ⫹ 1: 6 2

…16†

492

D.B. Parkinson / Advances in Engineering Software 30 (1999) 489–494

Fig. 2. Correlation functions and process sample for example 2.

Setting Rb(4a) ˆ 0 determines R 00b …0†: R 00b …0† ˆ

⫺10…1 ⫹ 4aR 0b …0††  : 251 113 0 ⫹ aR b 0 …0† a2 4 3

It will be assumed [19,23] that …17†

ai ˆ a1 pi⫺1 and

s2i ˆ s21 qi⫺1 for 3. The sum process As the simulated process is to be formed by the summation of N quadratic processes of the type described earlier, it is: S…t† ˆ

N X

bi …t†;

…18†

iˆ1

and its covariance is: CS …t† ˆ

N X

E‰bi …0†2 ŠRbi …t†

…20†

…19†

iˆ1

as cross-correlations of the bi are zero by virtue of the independence of the (wj)i, see Eq. (1).

0 ⬍ p; q ⱕ 1:

…21†

When, using Eqs. (12), (14), (17), (20) and (21):   251 113 0 i⫺1 0 R …0†p qi⫺1 ⫹ a 2 X 16s1 N 4 3 1 b 2 ; E…S…0† † ˆ 5 iˆ1 …1 ⫹ 4a1 R 0b 0 …0†pi⫺1 †…1 ⫹ 4a1 R 0b …0†pi⫺1 † …22† R 0bi 0 …0†

R 0b 0 …0†

ˆ and where it has been assumed that R 0bi …0† ˆ R 0b …0† for all i and their values will be determined to provide a best fit to the given correlation function. As the process to be simulated is standard normal E(S(0) 2) ˆ 1 and hence:   2 3⫺1=2 251 113 0 i⫺1 i⫺1 0 ⫹ a q R …0†p N 6 16 X 7 4 3 1 b 7 s1 ˆ 6 ; 4 5 0 i⫺1 0 i⫺1 …1 ⫹ 4a R …0†p †…1 ⫹ 4a R …0†p † 5 iˆ1

1

b0

1

b

…23†

D.B. Parkinson / Advances in Engineering Software 30 (1999) 489–494

493

Fig. 3. Correlation functions and process sample for example 3.

and the sum process correlation function is:   251 113 0 i⫺1 0 ⫹ a R …0†p qi⫺1 Rbi …t† N 16 2 X 4 3 1 b s ; RS …t† ˆ 5 1 iˆ1 …1 ⫹ 4a1 R 0b 0 …0†pi⫺1 †…1 ⫹ 4a1 R 0b …0†pi⫺1 † …24† while the crossing rate parameter g (see Refs. [15,18]) is given by: " #1=2 R 00s …0† gˆ ⫺ 2p 4s 1 ˆ a1

"

N 1 X qi⫺1 p iˆ1 …1 ⫹ 4a1 R 0b 0 …0†pi⫺1 †p2…i⫺1†

#1=2 :

…25†

4. Process simulation A set of points is selected from the positive time domain of the given correlation function and the function RS (Eq. (24)) is fitted to them by least squares optimisation [24] for a chosen value of the number of sub-processes, N. This then

determines the unknown parameters a1 ; p; q; R 0b 0 …0† and R 0b …0† which are the variables of optimisation and hence, s 1 (Eq. (23)), R 00b 0 …0† (Eq. (12)), R 00b …0† (Eq. (17)) and g (Eq. (25)). In the course of the optimisation procedure the sub-processes Rbi(t) are computed using the appropriate part of Eq. (16). With the parameters determined as mentioned earlier, each process bi(t) may be simulated using Eq. (1), with the wj determined from a standard normal random number generator and k from a uniform random number generator for the range zero to one. The resulting value of the sum process is then given by Eq. (18) at any time t and for the chosen value of N. As all the simulation parameters, except the random numbers are determined prior to the process simulation a fast simulation with minimal data storage requirements (three arrays of size N for the wj values) is achieved and it is possible to simulate selected regions or points of the process. Vector processes may be simulated by using the aforementioned as the component processes with specified cross correlations [19,22]. Non-normal processes may also be obtained using transformation methods together with the normal simulation described earlier [17,18].

494

D.B. Parkinson / Advances in Engineering Software 30 (1999) 489–494

Table 1 Process parameters Example number 1 2 3

A1 0.259 0.278 0.243

s1 0.015 0.013 0.026

p 0.899 0.935 0.953

5. Examples Three examples of correlation functions, fitted curves and corresponding process sample are provided in Figs. 1–3. The given function in the correlation graphs is shown by a continuous line, the approximating fitted curve by a dashed line, where only one curve is seen, the two curves are identical. Only the positive time halves of the functions are shown. The computations were carried out on a Sun/Unix system employing Uniras graphics. In each case 20 points, evenly spaced in time were selected from the given correlation function and used in the least squares curve fitting computation, using 30 subprocesses to form the sum process. The root mean square fitting errors in the three examples shown were 0.253, 0.227 and 0.032%, respectively, and the corresponding values of the main parameters a1, s 1, p, q, R 0b …0†; R 0b 0 …0† and g for each curve are listed in Table 1. In general it was found easier to fit a wider range of correlation function, using the second order processes described here, than for the original line process [18] and, in particular, an improved model for processes with correlation functions similar to example 3, was obtained. As explained in Ref. [18] the method may be used to obtain vector process with specified cross correlation and non-normal process by the use of transformation methods [17,18] and non-stationary processes, [2]. 6. Conclusions Previous work on stochastic process simulation using broken line processes has been extending to a wider range of correlation function and process type by the use of a quadratic sum process procedure. As storage of the whole process is unnecessary and the simulation procedure is fast, once the correlation function has been fitted, the method is particularly suited to uncertainty models (e.g. in reliability theory) requiring long duration simulations of single or vector processes associated with rare events. References [1] Borgman LE. Ocean wave simulation for engineering design. J Waterways Harbours Div. ASCE 1969;95(25):557–583. [2] Shinozuka M, Sato Y. Simulation of non-stationary random process. J. Engng Mech. Div. ASCE 1967;93(EMI):11–40.

R=b …0†

q 0.975 1.000 0.790

R=b 0 …0† ⫺13

⫺ 0.16 × 10 ⫺ 0.36 × 10 ⫺10 ⫺ 0.019

g

⫺ 0.17 × 10

⫺12

0 ⫺ 0.359

4.92 2.04 0.922

[3] Hon SN. Earthquake simulation models and their applications. Department of Civil Engineering, MIT, Rep Nr 68-17, 1968. [4] Davis BM, Hagan R, Borgman LE. A program for the FFT simulation of realisations from a one-dimensional random function with known covariance. Comput. Geosci. 1981;7:199–206. [5] Nishimura A. Stochastic consideration of the design live load of steel highway bridges. Mem. Faculty Engng Kobe Univ. 1957;4:37–49. [6] Mohapl J. A Stochastic model of water table elevation. Stochastic Hydrol. Hydrauls. 1998;12(4):223–245. [7] Arshem H. Stochastic optimisation of discrete event systems. Microelect. Reliab. 1996;36(10):1357–1368. [8] Ditlevsen O. Gaussian outcrossings from safe convex polyhedrons. J. Engng Mech. Div. ASCE 1983;109(1):127–148. [9] Madsen HO, Krenk S, Lund NC. Methods of structural safety. Englewood Cliffs, NJ: Prentice-Hall, 1986 pp. 147–194. [10] Rackwitz R. Failure rates for general systems including structural components. Reliab Engng 1984;9:229–242. [11] Parkinson DB. Multivariate uncertainty models with applications in design, reliability, maintainability and quality control. Doctoral Thesis, University of Liverpool, 1990. [12] Parkinson DB. Assessment and optimisation of dimensional tolerances. Computer Aided Design 1985;17(4):191–199. [13] Parkinson DB. Fast availability simulation. Reliab. Engng 1987;18(3):157–176. [14] Parkinson DB. Optimum sampling plans, based on post quality control reliability. Reliab. Engng System Safety 1988;21:59–75. [15] Morton-Firth CJ, Bray D. Predicting temporal fluctuations in an intracellular signalling pathway. J. Theoretical Biol. 1998;192(1):117– 128. [16] Wolters J. Stochastic dynamic properties of linear econometric models. Berlin: Springer, 1980. [17] Parkinson DB. First order reliability analysis employing translation systems. Engng Struct. 1978;1(1):31–40. [18] Grigoriu M. Crossings of non-Gaussian translation processes. J. Engng Mech. Div. ASCE 1984;110(4):610–620. [19] Ditlevsen O. Extremes and first passage times. Doctoral Thesis, Technical University of Denmark, Copenhagen, Denmark, 1971. [20] Christensen JK, Sorensen JD. Simulering af Gaussiske Processor na Datamaskine, Rep. 8007. Institute of Building Technology and Structural Engineering, AUC Aalborg, 1980, pp. 27–45. [21] Riply BD. Stochastic simulation. New York: Wiley, 1987, pp. 105– 110. [22] Parkinson DB. Simulation of stochastic processes with known correlation functions. Reliab. Engng System Safety 1993;40:213–220. [23] Hino M. Digital computer simulation of turbulent diffusion. International Association for Hydraulic Research, Congress Paper 2-30, Leningrad, 1965. [24] Greig DM. Optimisation. Harlow, UK: Longman, 1980, pp. 63–64. [25] Cramer H, Leadbetter MR. Stationary and related stochastic processes. New York: Wiley, 1967, 191–255. [26] Krenk S, Madsen PH. Stochastic response analysis. In: Thoft-Christensen P, editor. Reliability theory and its applications in structural and solid mechanics, 1983, pp. 103–172.