International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP) (12th : 2015)

Simulation of strongly non-Gaussian non-stationary stochastic processes utilizing Karhunen-Loeve expansion Kim, Hwanpyo; Shields, Michael D. 2015-07

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata


53032-Paper_139_Kim.pdf [ 491.65kB ]
JSON: 53032-1.0076040.json
JSON-LD: 53032-1.0076040-ld.json
RDF/XML (Pretty): 53032-1.0076040-rdf.xml
RDF/JSON: 53032-1.0076040-rdf.json
Turtle: 53032-1.0076040-turtle.txt
N-Triples: 53032-1.0076040-rdf-ntriples.txt
Original Record: 53032-1.0076040-source.json
Full Text

Full Text

12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015Simulation of Strongly Non-Gaussian Non-stationary StochasticProcesses utilizing Karhunen-Loeve ExpansionHwanpyo KimGraduate Student, Dept. of Civil Engineering, Johns Hopkins University, Baltimore, USAMichael D. ShieldsAssistant Professor, Dept. of Civil Engineering, Johns Hopkins University, Baltimore,USAABSTRACT: The simulation of non-stationary and non-Gaussian stochastic processes is a challengingproblem of considerable practical interest. Recently, Shields et al. have developed a class of concep-tually simple and efficient methods for simulation of non-Gaussian processes using translation processtheory (collectively referred to as the Iterative Translation Approximation Method - ITAM) that itera-tively upgrades the underlying Gaussian power spectral density function for simulation using the spectralrepresentation method. However, the currently existing ITAM method for generation of non-stationaryand non-Gaussian processes requires additional approximations in the estimation of the evolutionaryspectrum. An extension of the ITAM is proposed that utilizes the K-L expansion. The developed methoditeratively upgrades the covariance function directly and, in so doing, avoids the complex and non-uniqueinverse problem of estimating the evolutionary spectrum from the non-stationary autocorrelation. The ap-plication of the method for a strongly non-Gaussian and non-stationary process with a prescribed targetnon-Gaussian correlation function is demonstrated.1. INTRODUCTIONProbabilistic analysis of complex problems involv-ing strong non-linearities and uncertainties repre-sented by non-Gaussian stochastic processes/fieldstypically requires Monte Carlo simulation. MonteCarlo simulation, which is the most robust methodavailable for these problems, requires accurate sim-ulation of stochastic processes. Well establishedmethods are readily available for the simulation ofstationary Gaussian random processes. Unfortu-nately, these stationary and Gaussian processes areoften not representative of reality. Instead, realiza-tions of non-Gaussian and non-stationary stochasticprocess are necessary for accurate analysis. How-ever, simulation of these processes present numer-ous challenges and consequently, the developmentof methods for their simulation has lagged consid-erably behind the stationary and Gaussian models.The Karhunen-Loeve (K-L) expansion and theSpectral Representation Method (SRM) are the twomost commonly employed methods for generationof random processes in mechanics. The K-L ex-pansion (Ghanem and Spanos (1991), Huang et al.(2001)) affords an optimally succinct representa-tion of the stationary and non-stationary processesfrom a specified covariance function and its im-plementation is straightforward for Gaussian pro-cesses. Phoon et al. (2002) and Phoon et al.(2005) developed K-L-based methods for the sim-ulation of non-Gaussian and non-stationary pro-cesses that iteratively update the probability densityfunction’s of the K-L random variables. Further-more, Sakamoto and Ghanem (2002a,b) proposeda technique for simulating non-stationary and non-Gaussian stochastic processes that utilizes a poly-nomial chaos expansion (PCE) of the K-L randomvariables. Another class of algorithms for generat-ing non-Gaussian processes couples the SRM, (Shi-nozuka and Jan (1972); Shinozuka and Deodatis(1991)) with translation process theory, (Grigoriu(1995)). Several developments along these lineshave been proposed over the past 25 years includ-112th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015ing works by Yamazaki and Shinozuka (1988), De-odatis and Micaletti (2001), and Bocchini and De-odatis (2008) among others. Most recently, Shieldset al. (2011) and Shields and Deodatis (2013a)proposed the Iterative Translation ApproximationMethod (ITAM) for simulation of stationary non-Gaussian processes that upgrades the underlyingGaussian power spectral density function in a sim-ple and efficient manner. This method has beenextended to the simulation of non-stationary andnon-Gaussian processes by Shields and Deodatis(2013b).For strongly non-Gaussian and non-stationaryprocesses, each of the aforementioned methodspossesses notable drawbacks. K-L based simula-tion methods are subject to the Central Limit The-orem. Therefore, in the limit as the number of K-L random variables grows large the process tendstoward Gaussian regardless of the distribution ofthe K-L variables themselves. Moreover, even fora small number of K-L variables, the process of-ten cannot maintain the strongly non-Gaussian fea-tures that are desired. SRM-based methods, on theother hand, require significant approximations inorder to produce non-stationary and non-Gaussianprocesses (Shields and Deodatis (2013b)). Mostnotably, no unique inverse exists to compute theevolutionary spectrum (ES) as defined by Priest-ley (1965) from the non-stationary autocorrelationfunction. Thus, approximations must be introducedthat can cause significant errors for strongly non-Gaussian and non-stationary processes.In this work, an extension of the ITAM that uti-lizes the K-L expansion for simulation is proposed.The method developed in this study has the fol-lowing advantages for simulation of non-stationaryand non-Gaussian processes. First, estimation ofthe evolutionary spectrum is avoided; thus allevi-ating the associated approximations and allowingdirect simulation from the covariance matrix. Sec-ond, generated realizations are produced directlyfrom Gaussian K-L random variables with the re-sulting Gaussian process translated to match thetarget marginal non-Gaussian probability functionexactly. Lastly, It is conceptually simple, straight-forward to implement, and converges rapidly (usu-ally in 10 iterations or less). The applicationof the method for a strongly non-Gaussian andnon-stationary process with prescribed target non-Gaussian covariance is demonstrated.2. KARHUNEN-LOEVE EXPANSIONConsider a random process A(t,θ) defined onthe probability space (Ω,σ ,P) over the domains, t ∈ [a,b] with zero mean and covariance functionC(s, t) possessing finite variance. The process canbe expressed according to the K-L expansion asA(t,θ) =∞∑i=1√λiζi(θ) fi(t) (1)where λi and fi(x) are the eigenvalues and eigen-functions of the covariance C(s, t) determined bysolving the Fredholm integral equation of the sec-ond kind given by∫ baC(s, t) fi(s)ds = λi fi(t) (2)The parameters ζi are a set of uncorrelated ran-dom variables with zero mean E[ζi] = 0 and cor-relation E[ζiζ j] = δi j. The K-L expansion can beconveniently implemented for simulation by trun-cating a countable number, M, of its terms asA˜(t,θ) =M∑i=1√λiζi(θ) fi(x) (3)For Gaussian processes, the random variables ζifollow the standard Gaussian distribution. How-ever, the distribution of ζi for non-Gaussian pro-cesses can be difficult to determine and may requirecomplex dependence structure in general (Grigoriu(2010)).3. TRANSLATION PROCESS THEORYTranslation process theory (Grigoriu (1995)) repre-sents a non-Gaussian process Y (t) through the non-linear transformation of a Gaussian process, X(t),as follows:Y (t) = g(X(t)) (4)where g(·) is a general nonlinear transformation.The so-called standard translation is defined suchthat g(·) = F−1 · Φ(·) where F(·) and Φ(·) are212th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015the arbitrary non-Gaussian and standard Gaussianmarginal cumulative distribution functions, respec-tively.The extension of translation process theory fornon-stationary processes was proposed by Ferranteet al. (2005) where the mapping in Eq. (4) is gener-alized to incorporate a time-dependent transforma-tion asY (t) = g(X(t), t) = F−1 · {Φ(X(t)), t} (5)where X(t) is a non-stationary Gaussian processand Y (t) is a non-stationary and non-Gaussian pro-cess with time-varying marginal cumulative distri-bution function F(·, t).Utilizing non-stationary translation process the-ory, the correlation R(s, t) of the non-stationary andnon-Gaussian process can be computed from thestandard Gaussian correlation asR(s, t) = µ(s)µ(t)+√C(s,s)C(t, t)ξ (s, t)=∞∫−∞∞∫−∞g(u,s)g(v, t)φ{u,v;ρ(s, t)}dudv(6)where µ(t), C(t, t), and ξ (s, t) are the mean,variance, and correlation coefficient of the non-Gaussian process Y (t) and φ(s, t) denotes the jointGaussian PDF with correlation coefficient, ρ(s, t).A general challenge with translation process the-ory is that, given a prescribed covariance functionand associated non-Gaussian marginal probabilitydensity function, it is not always possible to iden-tify a corresponding underlying Gaussian correla-tion function. In other words, the inverse of Eq. (6)does not always exist. In such cases, the prescribedcovariance function and marginal probability den-sity function are said to be incompatible with thetranslation process model. However, given the con-venience of the translation model and the lack ofcomparable methods for modeling non-Gaussianprocesses, it is often desirable to identify an un-derlying Gaussian process that, when mapped tothe non-Gaussian distribution using Eq. (6), pro-duces a non-Gaussian correlation function that isas close as possible to the prescribed non-Gaussiancorrelation function. This is the motivation for thepreviously developed Iterative Translation Approx-imation Method (ITAM) described in the followingsection and its improvement in this work.4. ITERATIVE TRANSLATION APPROXI-MATION METHODThe ITAM, originally developed by Shields et al.(2011) and extended for non-Gaussian and non-stationary processes by Shields and Deodatis(2013b), utilizes a simple and efficient iterativemethod to upgrade the underlying Gaussian powerspectral density (stationary) or evolutionary spec-trum (non-stationary). In the non-stationary casethe ITAM upgrades the underlying Gaussian evolu-tionary spectrum asS(i+1)G (ω, t) =[STN(ω, t)S(i)N (ω, t)]βS(i)G (ω, t) (7)where S(i)G (ω, t) and S(i)N (ω, t) are the underlyingGaussian and the computed non-Gaussian ES at it-eration i, respectively and β is parameter utilized tooptimize the convergence.The challenge associated with this method isthat it requires estimation of the evolutionary spec-trum from the non-stationary correlation function.However, in general the evolutionary spectrum isnot uniquely defined for a given non-stationaryautocorrelation (Priestley (1965); Benowitz et al.(2014)). Thus, the ITAM proposed by Shields andDeodatis (2013b) requires an estimation that maycause significant errors in the case of strongly non-Gaussian and non-stationary processes. Recently,Benowitz et al. (2014) has investigated the unique-ness of this inversion and the results indicate thatit may be possible to obtain a unique evolutionaryspectrum under specific conditions, but the com-putational cost required to estimate this evolution-ary spectrum with reasonable accuracy is extremelyhigh.5. PROPOSED ITAM WITH K-L EXPAN-SIONA new ITAM utilizing the K-L expansion for sim-ulation is proposed in this work. This new ITAMiteratively upgrades the underlying non-stationaryautocorrelation function (ACF) directly (removing312th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015approximations associated with evolutionary spec-trum estimation). The following procedure de-scribes the proposed methodology for finding thecompatible underlying Gaussian ACF.1. Initialize underlying Gaussian ACF, R(0)G (s, t).2. Upgrade the underlying Gaussian ACF.3. Find the nearest positive semi-definite ACF tothe upgraded ACF.4. Compute the non-Gaussian ACF at iteration i,R(i)N (s, t), by non-stationary translation processtheory.5. Check the difference between the computedand target non-Gaussian ACF. If convergence,proceed to the next step. Otherwise, iterateback to step 2.6. Simulate the processes using the K-L expan-sion and translate process theory.When an incompatible pair of marginal non-Gaussian non-stationary PDF and ACF is pre-scribed, the initial underlying Gaussian ACF canbe defined arbitrarily - as long as it is a properlydefined autocorrelation function (i.e. is symmetricand positive semi-definite). In practice, it is use-ful to reduce calculation time by setting the initialunderlying Gaussian ACF equal to the target non-Gaussian ACF.The proposed upgrading method in step 2 is simi-lar to that of the SRM-based ITAM given in Eq. (7).Specifically, the underlying Gaussian normalizedACF, ρ(s, t) is upgraded asρ(i+1)(s, t) =[ξ T (s, t)ξ (i)(s, t)]ρ(i)(s, t) (8)where ξ (i)(s, t) is the computed non-Gaussian ACFat iteration i and ξ T (s, t) is the prescribed targetstandardized non-Gaussian ACF. However, thereare a few noteworthy differences with the SRM-based ITAM. Most notably, the evolutionary spec-trum is a strictly positive quantity while the ACF isnot. This has two significant consequences. First,the exponent β is removed in the proposed methodbecause iterations using a non-integer exponent willproduce imaginary numbers for negative ratios inthe parentheses. Second, the resulting ACF fromthe iterations in Eq. (8) is not necessary positivesemi-definite and is therefore not a valid autocor-relation function. To correct this, it is necessaryto identify the nearest positive semi-definite ACF(step 3). This is achieved using the method pro-posed by Qi and Sun (2006) who use a quadrati-cally convergent newton method whose quadraticconvergence has been proven.Finally, in step 5 the relative difference be-tween the computed and target non-Gaussian non-stationary ACF asε(i) = 100√√√√∑N−1n=0 ∑N−1m=0[ξ (i)(sn, tm)−ξ T (sn, tm)]2∑N−1n=0 ∑N−1m=0[ξ T (sn, tm)]2(9)When the value of the relative difference stabi-lizes to a constant value, the iterations are stoppedand the K-L expansion is used to generate theGaussian process from the converged underlyingGaussian ACF using Eq. (3). Then, the generatedGaussian sample function is translated to the non-Gaussian distribution using Eq. (5).6. NUMERICAL EXAMPLETo demonstrate the capabilities of the improvedITAM, consider the simulation of the strongly non-Gaussian and non-stationary process with targetnon-stationary covariance in Figure 1 is given byC(s, t) = min(s, t) · cos[4pi(s− t)] (10)Figure 1: Target non-stationary covariance.412th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015The process possesses a non-stationary "U-shaped" beta marginal distribution with cumulativedistribution function given byF(y; p,q) =Γ(p+q)Γ(p)Γ(q)∫ u0zp−1(1− z)q−1dz (11)where u = y−yminymax−ymin with upper and lower boundsymin = µb(t)+σb(t)√p(p+q+1)q and ymax = µb(t)−σb(t)√q(p+q+1)p , and Γ(·) is the gamma functionwith parameters p = 0.342 and q = 0.528. Themean µb(t) = 0 and the variance σ2b (t) =C(t, t) = t.The shape of this “U-shaped” beta PDF is shown inFigure 2 along with the standard normal pdf.x-6 -4 -2 0 2 4 6PDF00. Beta at t=0.25U Beta at t=0.5U Beta at t=0.75U Beta at t=1NormalFigure 2: Strongly non-Gaussian and non-stationary"U-shaped" beta marginal PDF shown with the stan-dard normal distributionThe converged non-Gaussian covariance at sev-eral time instants are shown in Figure 3. Despite thestrongly non-Gaussian and non-stationary nature ofthe target process, the computed non-Gaussian andnon-stationary covariance converges to the targetwith high accuracy for all values of t. The con-verged results are calculated in only 5 iterations.Across the entire time domain, the maximum rela-tive difference is 5.58% with the converged resultsmaintaining the shape and magnitude of the targetcovariance very well.s0 0.2 0.4 0.6 0.8 1Covariance-1-0.500.51 TargetComputed ϵ=2.866%(a) t = 1s0 0.2 0.4 0.6 0.8 1Covariance-1-0.500.51 TargetComputed ϵ=2.675%(b) t = 0.75s0 0.2 0.4 0.6 0.8 1Covariance-1-0.500.51 TargetComputed ϵ=3.309%(c) t = 0.5s0 0.2 0.4 0.6 0.8 1Covariance-1-0.500.51 TargetComputed ϵ=5.581%(d) t = 0.25Figure 3: Comparison of the computed non-Gaussianand non-stationary covariance with the target at differ-ent times7. COMPARISON WITH THE ORIGINALSRM-ITAMThe performance of the SRM-based ITAM and theproposed ITAM is demonstrated using the strongly512th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015non-stationary target ES given asS(t,ω) = 14√piexp[−(ω−ω0(t)2)2](12)where ω0(t) = 10+40t. Plots of the target ES andthe corresponding covariance are presented in Fig-ure 4. The prescribed PDF is the "U-shaped" betaPDF in Eq. (11) with constant mean µb(t) = 0 andvariance σ2b (t) = 1.(a) Target non-stationary evolutionary spectrum(b) Target non-stationary covarianceFigure 4: Target non-stationary evolutionary spectrumand covarianceNumerical results based on the two differentITAMs are presented on Table 1. As the originalTable 1: Relative difference of the computed non-Gaussian non-stationary covariance between the SRM-ITAM and KL-ITAMRelative difference (%)Time SRM-ITAM KL-TAMs = 0 17.81 1.086s = 1 12.11 5.317ITAM with SRM upgrades the underlying Gaus-sian ES, which is difficult to estimate for stronglynon-stationary cases and requires computation of apseudo-ES (Shields and Deodatis (2013a)), it haslarger relative difference throughout time. How-ever, the presented KL-ITAM, which upgrades theunderlying correlation function directly, approxi-mates the target covariance better than the previousmethod. The relative difference of the KL-basedITAM does not exceed 6%, the other difference isalmost 15% over the time. These differences willbe more significant for stronger non-Gaussian andnon-stationary targets.8. CONCLUSIONSAn improvement to the Iterative Translation Ap-proximation Method for generation of stronglynon-Gaussian and non-stationary process is pro-posed which utilizes direct upgrading on the non-stationary autocorrelation function and simulationwith the K-L expansion. The method improvesupon existing K-L based simulation methods bypairing them with translation process theory in or-der to match marginal non-Gaussian probabilitydensity functions exactly. It also improves uponthe existing SRM-based ITAM simulation methodsby removing the need to estimate the evolutionaryspectrum which can produce significant errors. Themethod is shown to be highly accurate, straightfor-ward to implement, and converge rapidly - even forprocesses that are very difficult to simulate.9. REFERENCESBenowitz, B., Shields, M., and Deodatis, G. (2014).“Determining evolutionary spectra from non-stationary autocorrelation functions.” ProbabilisticEngineering Mechanics.612th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015t0 0.5 1 1.5 2Covariance-1-0.500.51TargetSRM-ITAMKL-ITAM(a) s = 0t0 0.5 1 1.5 2Covariance-1-0.500.51 TargetSRM-ITAMKL-ITAM(b) s = 1Figure 5: Comparison of the computed non-Gaussianand non-stationary covariance with the target betweenSRM-ITAM and KL-ITAMBocchini, P. and Deodatis, G. (2008). “Critical re-view and latest developments of a class of simu-lation algorithms for strongly non-gaussian randomfields.” Probabilistic Engineering Mechanics, 23(4),393–407.Deodatis, G. and Micaletti, R. C. (2001). “Simulationof highly skewed non-gaussian stochastic processes.”Journal of engineering mechanics, 127(12), 1284–1295.Ferrante, F., Arwade, S., and Graham-Brady, L. (2005).“A translation model for non-stationary, non-gaussianrandom processes.” Probabilistic engineering me-chanics, 20(3), 215–228.Ghanem, R. G. and Spanos, P. D. (1991). Stochasticfinite elements: a spectral approach, Vol. 387974563.Springer.Grigoriu, M. (1995). Applied non-Gaussian processes:Examples, theory, simulation, linear random vibra-tion, and MATLAB solutions. PTR Prentice Hall Up-per Saddle River, NJ.Grigoriu, M. (2010). “Probabilistic models for stochas-tic elliptic partial differential equations.” Journal ofComputational Physics, 229(22), 8406–8429.Huang, S., Quek, S., and Phoon, K. (2001). “Conver-gence study of the truncated karhunen–loeve expan-sion for simulation of stochastic processes.” Interna-tional journal for numerical methods in engineering,52(9), 1029–1043.Phoon, K., Huang, H., and Quek, S. (2005). “Simulationof strongly non-gaussian processes using karhunen–loève expansion.” Probabilistic engineering mechan-ics, 20(2), 188–198.Phoon, K., Huang, S., and Quek, S. (2002). “Simula-tion of second-order processes using karhunen–loeveexpansion.” Computers & structures, 80(12), 1049–1060.Priestley, M. B. (1965). “Evolutionary spectra and non-stationary processes.” Journal of the Royal StatisticalSociety. Series B (Methodological), 204–237.Qi, H. and Sun, D. (2006). “A quadratically convergentnewton method for computing the nearest correlationmatrix.” SIAM journal on matrix analysis and appli-cations, 28(2), 360–385.Sakamoto, S. and Ghanem, R. (2002a). “Polyno-mial chaos decomposition for the simulation of non-gaussian nonstationary stochastic processes.” Journalof engineering mechanics, 128(2), 190–201.Sakamoto, S. and Ghanem, R. (2002b). “Simulation ofmulti-dimensional non-gaussian non-stationary ran-dom fields.” Probabilistic Engineering Mechanics,17(2), 167–176.Shields, M. and Deodatis, G. (2013a). “Estimation ofevolutionary spectra for simulation of non-stationaryand non-gaussian stochastic processes.” Computers &Structures, 126, 149–163.Shields, M. and Deodatis, G. (2013b). “A simple andefficient methodology to approximate a general non-gaussian stationary stochastic vector process by atranslation process with applications in wind veloc-ity simulation.” Probabilistic Engineering Mechan-ics, 31, 19–29.Shields, M., Deodatis, G., and Bocchini, P. (2011). “Asimple and efficient methodology to approximate a712th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015general non-gaussian stationary stochastic process bya translation process.” Probabilistic Engineering Me-chanics, 26(4), 511–519.Shinozuka, M. and Deodatis, G. (1991). “Simulation ofstochastic processes by spectral representation.” Ap-plied Mechanics Reviews, 44(4), 191–204.Shinozuka, M. and Jan, C.-M. (1972). “Digital simula-tion of random processes and its applications.” Jour-nal of sound and vibration, 25(1), 111–128.Yamazaki, F. and Shinozuka, M. (1988). “Digital gen-eration of non-gaussian stochastic fields.” Journal ofEngineering Mechanics, 114(7), 1183–1197.8


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items