git clone https://git.igankevich.com/waves-16-arma.git
Log | Files | Refs

commit ac0a4ffbef95fd97ca086717d2d06bed4ee1a313
parent e2a1c2f1180ac53c44caf6a4911ded3e0c9a8d7d
Author: Ivan Gankevich <igankevich@ya.ru>
Date:   Sun, 28 May 2017 09:44:16 +0300

Add related work and governing equations.

arma.org | 396++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
preamble.tex | 30++++++++++++++++++++++++++++--
refs.bib | 82+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--
3 files changed, 501 insertions(+), 7 deletions(-)

diff --git a/arma.org b/arma.org @@ -4,7 +4,7 @@ #+LATEX_CLASS: scrartcl #+LATEX_CLASS_OPTIONS: #+LATEX_HEADER_EXTRA: \input{preamble} -#+OPTIONS: H:2 num:0 todo:nil toc:nil +#+OPTIONS: H:5 num:0 todo:nil toc:nil #+begin_abstract Simulation of sea waves is a problem appearing in the framework of developing @@ -73,7 +73,7 @@ isotropy within the marine environment. In this case eq.\nbsp{}eqref:eq-continui reduces to \begin{equation} \label{eq-continuity-2} - \mathoperator{div}\vec{V} = + \text{div}\vec{V} = \frac{\partial{V_x}}{\partial{x}} + \frac{\partial{V_y}}{\partial{y}} + \frac{\partial{V_z}}{\partial{z}} = 0 @@ -145,13 +145,403 @@ becomes an even more significant issue in a nonlinear computation where the wave model is even more complex. Thus identifying a significantly less time intensive method for modelling the ambient ocean-wave environment has the potential for significantly speeding the total simulation process. - * Related work +In\nbsp{}cite:spanos1982arma ARMA model is used to generate time series spectrum of +which is compatible with Pierson---Moskowitz (PM) approximation of ocean wave +spectrum. The authors carry out experiments for one-dimensional AR, MA and ARMA +models. They mention excellent agreement between target and initial spectra and +higher performance of ARMA model compared to models based on summing large +number of harmonic components with random phases. The also mention that in order +to reach agreement between target and initial spectrum MA model require lesser +number of coefficients than AR model. In\nbsp{}cite:spanos1996efficient the authors +generalise ARMA model coefficients determination formulae for multi-variate +(vector) case. + +One thing that distinguishes present work with respect to afore-mentioned ones +is the study of three-dimensional (2D in space and 1D in time) ARMA model, which +is mostly a different problem. +1. Yule---Walker system of equations, which are used to determine AR + coefficients, has complex block-block structure. +2. Optimal model order (in a sense that target spectrum agrees with initial) is + determined manually. +3. Instead of PM spectrum, analytic formulae for standing and propagating + waves ACF are used as the model input. +4. Three-dimensional wavy surface should be compatible with real ocean surface + not only in terms of spectral characteristics, but also in the shape of wave + profiles. So, model verification includes distributions of various parameters + of generated waves (lengths, heights, periods etc.). +Multi-dimensionality of investigated model not only complexifies the task, but +also allows to carry out visual validation of generated wavy surface. It is the +opportunity to visualise output of the programme that allowed to ensure that +generated surface is compatible with real ocean surface, and is not abstract +multi-dimensional stochastic process that is real only statistically. + * Three-dimensional ARMA process as a sea wave simulation model + +Another approach to modelling wind waves is possible in terms of the +representation of a stochastic moving surface as a linear transformation of +white noise with memory. These methods are one of the most popular ways of +modelling stationary ergodic Gaussian random processes with given correlation +characteristics (Box, et al., 2008). However, these methods have were not used +to simulate wind waves for a long time. The first attempts to model +two-dimensional disturbances were undertaken in the early 70's (cf. Kostecki, +1972), and the impetus for this was the development of the resonance theory of +waves in wind. However, the formal mathematical framework was developed by +Gurgenidze & Trapeznikov (1988) and Rozhkov & Trapeznikov (1990). They built a +one-dimensional model of ocean waves \(\zeta(t)\), on the basis of an +autoregressive-moving average (ARMA) model: +\begin{equation} + \zeta_{\vec i} + = + \sum\limits_{\vec j = \vec 0}^{\vec N} + \Phi_{\vec j} \zeta_{\vec i - \vec j} + + + \sum\limits_{\vec j = \vec 0}^{\vec M} + \Theta_{\vec j} \epsilon_{\vec i - \vec j} + , + \label{eq-arma-process} +\end{equation} +where \(\zeta\)\nbsp{}--- wave elevation, \(\Phi\)\nbsp{}--- AR process +coefficients, \(\Theta\)\nbsp{}--- MA process coefficients, +\(\epsilon\)\nbsp{}--- white noise with Gaussian distribution, +\(\vec{N}\)\nbsp{}--- AR process order, \(\vec{M}\)\nbsp{}--- MA process order, +and \(\Phi_{\vec{0}}\equiv{0}\), \(\Theta_{\vec{0}}\equiv{0}\). Here arrows +denote multi-component indices with a component for each dimension. In general, +any scalar quantity can be a component (temperature, salinity, concentration of +some substance in water etc.). Equation parameters are AR and MA process +coefficients and order. + +Any ARMA process can be uniquely represented as a process moving average and +autoregression process of general infinite order (Gurgenidze & Trapeznikov, +1988), and the parameters of the spectral representation are defined by the rule +of division of power series (in a rational factorized form, Rozhkov & +Trapeznikov, 1990): +\begin{equation*} + S(\omega) = + \frac{\Delta\sigma^2}{\pi} + \frac{\prod\limits_m (1 - z_m e^{-im\omega\delta})(1 - z_m e^{im\omega\delta})} + {\prod\limits_n (1 - p_n e^{-in\omega\delta})(1 - p_n e^{in\omega\delta})}, +\end{equation*} +where \(z_m\) and \(p_n\) are the zeros of numerator (MA), and denominator (AR), +respectively, which form a pair of mutually conjugate numbers. If some of the +zeros are located near the unit circle, then the spectral density will have +pronounced dips. + ** Autoregressive (AR) process +AR process is ARMA process with only one random impulse instead of theirs +weighted sum: +\begin{equation} + \zeta_{\vec i} + = + \sum\limits_{\vec j = \vec 0}^{\vec N} + \Phi_{\vec j} \zeta_{\vec i - \vec j} + + + \epsilon_{i,j,k} + . + \label{eq-ar-process} +\end{equation} +The coefficients \(\Phi\) are calculated from ACF via three-dimensional +Yule---Walker equations, which are obtained after multiplying both parts of the +previous equation by \(\zeta_{\vec{i}-\vec{k}}\) and computing the expected value. +Generic form of YW equations is +\begin{equation} + \label{eq-yule-walker} + \gamma_{\vec k} + = + \sum\limits_{\vec j = \vec 0}^{\vec N} + \Phi_{\vec j} + \text{ }\gamma_{\vec{k}-\vec{j}} + + + \Var{\epsilon} \delta_{\vec{k}}, + \qquad + \delta_{\vec{k}} = + \begin{cases} + 1, \quad \text{if } \vec{k}=0 \\ + 0, \quad \text{if } \vec{k}\neq0, + \end{cases} +\end{equation} +where \(\gamma\)\nbsp{}--- ACF of process \(\zeta\), \(\Var{\epsilon}\)\nbsp{}--- white noise +variance. Matrix form of three-dimensional YW equations, which is used in the +present work, is +\begin{equation*} + \Gamma + \left[ + \begin{array}{l} + \Phi_{\vec 0}\\ + \Phi_{0,0,1}\\ + \vdotswithin{\Phi_{\vec 0}}\\ + \Phi_{\vec N} + \end{array} + \right] + = + \left[ + \begin{array}{l} + \gamma_{0,0,0}-\Var{\epsilon}\\ + \gamma_{0,0,1}\\ + \vdotswithin{\gamma_{\vec 0}}\\ + \gamma_{\vec N} + \end{array} + \right], + \qquad + \Gamma= + \left[ + \begin{array}{llll} + \Gamma_0 & \Gamma_1 & \cdots & \Gamma_{N_1} \\ + \Gamma_1 & \Gamma_0 & \ddots & \vdotswithin{\Gamma_0} \\ + \vdotswithin{\Gamma_0} & \ddots & \ddots & \Gamma_1 \\ + \Gamma_{N_1} & \cdots & \Gamma_1 & \Gamma_0 + \end{array} + \right], +\end{equation*} +where \(\vec N = \left( p_1, p_2, p_3 \right)\) and +\begin{equation*} + \Gamma_i = + \left[ + \begin{array}{llll} + \Gamma^0_i & \Gamma^1_i & \cdots & \Gamma^{N_2}_i \\ + \Gamma^1_i & \Gamma^0_i & \ddots & \vdotswithin{\Gamma^0_i} \\ + \vdotswithin{\Gamma^0_i} & \ddots & \ddots & \Gamma^1_i \\ + \Gamma^{N_2}_i & \cdots & \Gamma^1_i & \Gamma^0_i + \end{array} + \right] + \qquad + \Gamma_i^j= + \left[ + \begin{array}{llll} + \gamma_{i,j,0} & \gamma_{i,j,1} & \cdots & \gamma_{i,j,N_3} \\ + \gamma_{i,j,1} & \gamma_{i,j,0} & \ddots &x \vdotswithin{\gamma_{i,j,0}} \\ + \vdotswithin{\gamma_{i,j,0}} & \ddots & \ddots & \gamma_{i,j,1} \\ + \gamma_{i,j,N_3} & \cdots & \gamma_{i,j,1} & \gamma_{i,j,0} + \end{array} + \right], +\end{equation*} +Since \(\Phi_{\vec 0}\equiv0\), the first row and column of \(\Gamma\) can be +eliminated. Matrix \(\Gamma\) is block-toeplitz, positive definite and symmetric, +hence the system is efficiently solved by Cholesky decomposition, which is +particularly suitable for these types of matrices. + +After solving this system of equations white noise variance is estimated from +eqref:eq-yule-walker by plugging \(\vec k = \vec 0\): +\begin{equation*} + \Var{\epsilon} = + \Var{\zeta} + - + \sum\limits_{\vec j = \vec 0}^{\vec N} + \Phi_{\vec j} + \text{ }\gamma_{\vec{j}}. +\end{equation*} ** Moving average (MA) process +MA process is ARMA process with \(\Phi\equiv0\): +\begin{equation} + \zeta_{\vec i} + = + \sum\limits_{\vec j = \vec 0}^{\vec M} + \Theta_{\vec j} \epsilon_{\vec i - \vec j} + . + \label{eq-ma-process} +\end{equation} +MA coefficients \(\Theta\) are defined implicitly via the following non-linear +system of equations: +\begin{equation*} + \gamma_{\vec i} = + \left[ + \displaystyle + \sum\limits_{\vec j = \vec i}^{\vec M} + \Theta_{\vec j}\Theta_{\vec j - \vec i} + \right] + \Var{\epsilon}. +\end{equation*} +The system is solved numerically by fixed-point iteration method via the +following formulae +\begin{equation*} + \Theta_{\vec i} = + -\frac{\gamma_{\vec 0}}{\Var{\epsilon}} + + + \sum\limits_{\vec j = \vec i}^{\vec M} + \Theta_{\vec j} \Theta_{\vec j - \vec i}. +\end{equation*} +Here coefficients \(\Theta\) are calculated from back to front: from +\(\vec{i}=\vec{M}\) to \(\vec{i}=\vec{0}\). White noise variance is estimated by +\begin{equation*} + \Var{\epsilon} = \frac{\gamma_{\vec 0}}{ + 1 + + + \sum\limits_{\vec j = \vec 0}^{\vec M} + \Theta_{\vec j}^2 + }. +\end{equation*} +Authors of\nbsp{}cite:box1976time suggest using Newton---Raphson method to solve this +equation with higher precision, however, this method does not work in three +dimensions. Using slower method does not have dramatic effect on the overall +programme performance, because the number of coefficients is small and most of +the time is spent generating wavy surface. ** Mixed autoregressive moving average (ARMA) process +:PROPERTIES: +:CUSTOM_ID: sec:how-to-mix-ARMA +:END: +Generally speaking, ARMA process is obtained by plugging MA generated wavy +surface as random impulse to AR process, however, in order to get the process +with desired ACF one should re-compute AR coefficients before plugging. There +are several approaches to "mix" AR and MA processes. +- The approach proposed in\nbsp{}cite:box1976time which involves dividing ACF into MA + and AR part along each dimension is not applicable here, because in three + dimensions such division is not possible: there always be parts of the ACF + that are not taken into account by AR and MA process. +- The alternative approach is to use the same (undivided) ACF for both AR and MA + processes but use different process order, however, then realisation + characteristics (mean, variance etc.) become skewed: these are characteristics + of the two overlapped processes. +For the first approach there is a formula to re-compute ACF for AR process, but +there is no such formula for the second approach. So, the best solution for now +is to simply use AR and MA process exclusively. + +** Process selection criteria for different wave profiles +:PROPERTIES: +:CUSTOM_ID: sec-process-selection +:END: +One problem of ARMA model application to ocean wave generation is that for +different types of wave profiles different processes /must/ be used: standing +waves are modelled by AR process, and propagating waves by MA process. This +statement comes from practice: if one tries to use the processes the other way +round, the resulting realisation either diverges or does not correspond to real +ocean waves. (The latter happens for non-invertible MA process, as it is always +stationary.) So, the best way to apply ARMA model to ocean wave generation is to +use AR process for standing waves and MA process for progressive waves. + +The other problem is inability to automatically determine optimal number of +coefficients for three-dimensional AR and MA processes. For one-dimensional +processes this can be achieved via iterative methods\nbsp{}cite:box1976time, but they +diverge in three-dimensional case. + +The final problem, which is discussed in [[#sec:how-to-mix-ARMA]], is inability to +"mix" AR and MA process in three dimensions. + +In practice some statements made for AR and MA processes in\nbsp{}cite:box1976time +should be flipped for three-dimensional case. For example, the authors say that +ACF of MA process cuts at \(q\) and ACF of AR process decays to nought infinitely, +but in practice making ACF of 3-dimensional MA process not decay results in it +being non-invertible and producing realisation that does not look like real +ocean waves, whereas doing the same for ACF of AR process results in stationary +process and adequate realisation. Also, the authors say that one +should allocate the first \(q\) points of ACF to MA process (as it often needed to +describe the peaks in ACF) and leave the rest points to AR process, but in +practice in case of ACF of a propagating wave AR process is stationary only for +the first time slice of the ACF, and the rest is left to MA process. + +To summarise, the only established scenario of applying ARMA model to ocean wave +generation is to use AR process for standing waves and MA process for +propagating waves. With new formulae for 3 dimensions a single mixed ARMA +process might increase model precision, which is one of the objectives of the +future research. + ** The shape of ACF for different types of waves +**** Analytic method of finding the ACF. +The straightforward way to find ACF for a given ocean wave profile is to apply +Wiener---Khinchin theorem. According to this theorem the autocorrelation \(K\) of +a function \(\zeta\) is given by the Fourier transform of the absolute square of +the function: +\begin{equation} + K(t) = \Fourier{\left| \zeta(t) \right|^2}. + \label{eq-wiener-khinchin} +\end{equation} +When \(\zeta\) is replaced with actual wave profile, this formula gives you +analytic formula for the corresponding ACF. + +For three-dimensional wave profile (2D in space and 1D in time) analytic formula +is a polynomial of high order and is best obtained via symbolic computation +programme. Then for practical usage it can be approximated by superposition of +exponentially decaying cosines (which is how ACF of a stationary ARMA process +looks like\nbsp{}cite:box1976time). + +**** Empirical method of finding the ACF. +However, for three-dimensional case there exists simpler empirical method which +does not require sophisticated software to determine shape of the ACF. It is +known that ACF represented by exponentially decaying cosines satisfies first +order Stokes' equations for gravity waves\nbsp{}cite:boccotti1983wind. So, if the +shape of the wave profile is the only concern in the simulation, then one can +simply multiply it by a decaying exponent to get appropriate ACF. This ACF does +not reflect other wave profile parameters, such as wave height and period, but +opens possibility to simulate waves of a particular non-analytic shape by +"drawing" their profile, then multiplying it by an exponent and using the +resulting function as ACF. So, this empirical method is imprecise but offers +simpler alternative to Wiener---Khinchin theorem approach; it is mainly useful +to test ARMA model. + +**** Standing wave ACF. +For three-dimensional plain standing wave the profile is given by +\begin{equation} + \zeta(t, x, y) = A \sin (k_x x + k_y y) \sin (\sigma t). + \label{eq-standing-wave} +\end{equation} +Find ACF via analytic method. Multiplying the formula by a decaying exponent +(because Fourier transform is defined for a function \(f\) that +\(f\underset{x\rightarrow\pm\infty}{\longrightarrow}0\)) yields +\begin{equation} + \zeta(t, x, y) = + A + \exp\left[-\alpha (|t|+|x|+|y|) \right] + \sin (k_x x + k_y y) \sin (\sigma t). + \label{eq-decaying-standing-wave} +\end{equation} +Then, apply 3D Fourier transform to both sides of the equation via symbolic +computation programme, fit the resulting polynomial to the following +approximation: +\begin{equation} + K(t,x,y) = + \gamma + \exp\left[-\alpha (|t|+|x|+|y|) \right] + \cos \beta t + \cos \left[ \beta x + \beta y \right]. + \label{eq-standing-wave-acf} +\end{equation} +So, after applying Wiener---Khinchin theorem we get initial formula but with +cosines instead of sines. This difference is important because the value of ACF +at \((0,0,0)\) equals to the ARMA process variance, and if one used sines the +value would be wrong. + +If one tries to replicate the same formula via empirical method, the usual way +is to adapt eqref:eq-decaying-standing-wave to match eqref:eq-standing-wave-acf. +This can be done either by changing the phase of the sine, or by substituting +sine with cosine to move the maximum of the function to the origin of +coordinates. + +**** Propagating wave ACF. +Three-dimensional profile of plain propagating wave is given by +\begin{equation} + \zeta(t, x, y) = A \cos (\sigma t + k_x x + k_y y). + \label{eq-propagating-wave} +\end{equation} +For the analytic method repeating steps from the previous two paragraphs yields +\begin{equation} + K(t,x,y) = + \gamma + \exp\left[-\alpha (|t|+|x|+|y|) \right] + \cos\left[\beta (t+x+y) \right]. + \label{eq-propagating-wave-acf} +\end{equation} +For the empirical method the wave profile is simply multiplied by a decaying +exponent without need to adapt the maximum value of ACF (as it is required for +standing wave). + +**** Comparison of studied methods. +To summarise, the analytic method of finding ocean wave's ACF reduces to the +following steps. +- Make wave profile decay when approaching \(\pm\infty\) by multiplying it by + a decaying exponent. +- Apply Fourier transform to the absolute square of the resulting equation using + symbolic computation programme. +- Fit the resulting polynomial to the appropriate ACF approximation. + +Two examples in this section showed that in case of standing and propagating +waves their decaying profiles resemble the corresponding ACFs with the exception +that the ACF's maximum should be moved to the origin to preserve simulated +process variance. Empirical method of finding ACF reduces to the following +steps. +- Make wave profile decay when approaching \(\pm\infty\) by multiplying it by + a decaying exponent. +- Move maximum value of the resulting function to the origin by using + trigonometric identities to shift the phase. + ** Evaluation ** Discussion * Determining wave pressures for discretely given wavy surface diff --git a/preamble.tex b/preamble.tex @@ -1,2 +1,28 @@ \usepackage{amsmath} -\usepackage{cite}- \ No newline at end of file +\usepackage{cite} +\usepackage{latexsym} % \Box macro +\usepackage{mathtools} % fancy dots in matrices + +% custom mathematical expressions +\newcommand{\Var}[1]{\sigma_{#1}^2} +\newcommand{\Fourier}[1]{\mathcal{F}\left\{#1\right\}} +\newcommand{\InverseFourier}[1]{\mathcal{F}^{-1}\left\{#1\right\}} +\newcommand{\Fun}[1]{\mathcal{D}_1\left(x,#1\right)} +\newcommand{\FunSecond}[1]{\mathcal{D}_2\left(x,#1\right)} +\newcommand{\FunThird}[1]{\mathcal{D}_2\left(x,#1\right)} +\newcommand{\FunThreeD}[1]{\mathcal{D}_3\left(x,y,#1\right)} +\newcommand{\Sinh}[1]{\cosh\left(#1\right)} +\newcommand{\SinhX}[1]{\sinh\left(#1\right)} + +\newcommand{\FourierY}[2]{\mathcal{F}_{#2}\!\left\{#1\right\}} +\newcommand{\InverseFourierY}[2]{\mathcal{F}^{-1}_{#2}\!\left\{#1\right\}} + +\newcommand{\FourierX}[3]{\mathcal{F}_{#2}\!\left\{#1\right\}\!\left(#3\right)} +\newcommand{\InverseFourierX}[3]{\mathcal{F}^{-1}_{#2}\!\left\{#1\right\}\!\left(#3\right)} + +% properly aligned version of sqrt for \zeta_y^2 +\newcommand{\SqrtZeta}[1]{\sqrt{\vphantom{\zeta_x^2}\smash[b]{#1}}} + +% wave vector +\newcommand{\Kvec}{\vec{k}} +\newcommand{\Kveclen}{\lvert\smash[b]{\Kvec}\rvert}+ \ No newline at end of file diff --git a/refs.bib b/refs.bib @@ -34,17 +34,94 @@ @book{kochin1964theoretical, title={Theoretical hydromechanics}, - author={Kochin, Nikola{\u\i}} and Iliia, A Kibel and Roze, Nikola{\u\i}}}, + author={Kochin, Nikolai and Iliia, A Kibel and Roze, Nikolai}, year={1964}, publisher={Interscience} } @article{beck2001modern, title={Modern computational methods for ships in a seaway}, - author={BECK, Robert F and REED, Arthur M and SCLAVOUNOS, Paul D and HUTCHISON, Bruce L}, + author={Beck, Robert F and Reed, Arthur M and Sclavounos, Paul D and Hutchison, Bruce L}, journal={Transactions-Society of Naval Architects and Marine Engineers}, volume={109}, pages={1--51}, year={2001}, publisher={Society of Naval Architects and Marine Engineers} } + +@article{spanos1983arma, + title={ARMA algorithms for ocean wave modeling}, + author={Spanos, PT D}, + journal={Journal of Energy Resources Technology}, + volume={105}, + number={3}, + pages={300--309}, + year={1983}, + publisher={American Society of Mechanical Engineers} +} + +@book{spanos1982arma, + title={ARMA Algorithms for Ocean Spectral Analysis}, + author={Spanos, Pol D}, + year={1982}, + publisher={University of Texas at Austin. Engineering Mechanics Research Laboratory} +} + +@article{mignolet1992simulation, + title={Simulation of homogeneous two-dimensional random fields: Part I—AR and ARMA models}, + author={Mignolet, Marc P and Spanos, Pol D}, + journal={Journal of applied mechanics}, + volume={59}, + number={2S}, + pages={S260--S269}, + year={1992}, + publisher={American Society of Mechanical Engineers} +} + +@article{spanos1992simulation, + title={Simulation of homogeneous two-dimensional random fields: Part II—MA and ARMA Models}, + author={Spanos, Pol D and Mignolet, Marc P}, + journal={Journal of applied mechanics}, + volume={59}, + number={2S}, + pages={S270--S277}, + year={1992}, + publisher={American Society of Mechanical Engineers} +} + +@article{spanos1996efficient, + title={Efficient iterative ARMA approximation of multivariate random processes for structural dynamics applications}, + author={Spanos, Pol D and Zeldin, BA}, + journal={Earthquake engineering \& structural dynamics}, + volume={25}, + number={5}, + pages={497--507}, + year={1996}, + publisher={Wiley Online Library} +} + +@phdthesis{zeldin1996representation, + title={Representation and synthesis of random fields: ARMA, Galerkin, and wavelet procedures}, + author={Zeldin, Boris A}, + year={1996}, + school={Rice University} +} + +@Book{ box1976time, + title = {Time series analysis: forecasting and control, revised + ed}, + author = {Box, George EP and Jenkins, Gwilym M}, + year = {1976}, + publisher = {Holden-Day} +} + +@Article{ boccotti1983wind, + title = {On wind wave kinematics}, + author = {Boccotti, Paolo}, + journal = {Meccanica}, + volume = {18}, + number = {4}, + pages = {205--216}, + year = {1983}, + publisher = {Springer} +}+ \ No newline at end of file