

\section{Introduction}

The last decade has seen a dramatic increase in the use of Monte Carlo methods for solving stochastic engineering problems. There are primarily two reasons for this increase. First, the computational power available today, even for a laptop computer, is formidable and steadily increasing. Second, the versatility of Monte Carlo methods make them very attractive as a way of obtaining solutions to stochastic problems. The drawback of Monte Carlo methods for a range of problems has been that the required numerical calculations may take days, weeks or even months to do. But this situation is changing, some numerical problems that required several days of computer time for their solution just a few years ago can now be solved in minutes or hours. This has really opened the door for the use of Monte Carlo-based methods for solving a wide array of stochastic engineering problems. In this chapter the focus is on adapting Monte Carlo methods for estimation of extreme values of stochastic processes encountered in various engineering disciplines.        




\section{Simulation of Stationary Stochastic Processes} \label{sec:sssp} \index{stochastic process!simulation}


The approach to the simulation of stationary stochastic processes
favored in this book, is the spectral representation method \citep{shde91}. The
main reason for this is its simplicity and transparency for
practical applications.

The procedure is described in Subsection~\ref{ssec:variancespectrum},
where Eq.~(\ref{eq:4.8}) is a key result. Useful background
information for this is contained in Subsection~\ref{ssec:realizations}. 
Combining this with Example~\ref{ex:51} provides a
hands-on guide of how to simulate realizations of a stationary
Gaussian process. 

The use of the fast Fourier transform (FFT) technique can
substantially speed up the production of realizations of
stationary Gaussian processes, cf. \cite{shino74,newl91}.




\subsection{Realizations of stochastic processes }\label{ssec:realizations} \index{stochastic process!realization}

Monte Carlo methods can be used for any input-output system
subjected to a stochastic process such that for every realization
of the input stochastic process a corresponding realization of the
output process can be calculated. The basic idea underlying the
Monte Carlo method is that of producing a sample of
output/response time histories from a sample of input/loading time
histories. This makes it possible to estimate various statistics
of the output/response process based on the available sample of
realizations. Hence, a key element in the Monte Carlo method is
therefore the realizations of a stochastic process. For their simplicity, the focus here is on a couple of widely used ways of representing stationary stochastic processes, and a discussion is given of how
realizations of such processes are generated, see also \cite{shja72,shde91}. 

A simple, but useful representation of a stationary stochastic process $X(t)$ can be obtained as follows,
\begin{equation}\label{eq:8.1}
X(t) = \sum_{j=1}^N \, \big\{ A_j \, \cos \omega_j t + B_j \, \sin
\omega_j t \big\} \, , \quad \- \infty < t < \infty, 
\end{equation}
where, for $j = 1, \ldots, N$, the $\omega_j$ are positive constants, the $A_j$ and the $B_j$ are random variables with
$\EE[A_j] = \EE[B_j] = 0$, $\EE[A_j^2] = \EE[B_j^2] = \sigma_j^2$, and
$\EE[A_j B_j] = 0$ for $j = 1, \ldots, N$. Also, $\EE[A_j A_k] =
\EE[B_j B_k] = \EE[A_j B_k] = 0$ for $j \ne k$. When these
conditions are satisfied, it can be shown 
that $X(t)$ is a (weakly) stationary process with,
\begin{equation}\label{eq:8.3}
m_X = \EE[X(t)] = 0, 
\end{equation}
and
\begin{equation}\label{eq:8.4}
C_X(\tau) = \EE[(X(t)-m_X)( X(t + \tau) - m_X)] = \sum_{j=1}^N \sigma_j^2 \, \cos \omega_j
\tau. 
\end{equation}

The first step of the procedure is to specify
what kind of random variables $A_j$ and $B_j$ are assumed to be. A
common choice in many cases is to assume that these variables are
independent and normally distributed with zero mean value and
known standard deviations $\sigma_j$. A realization of $X(t)$ is
then obtained when a set of outcomes of the
random variables $A_j$ and $B_j$ has been generated. In this connection, it is useful
to make the following observation: if $\tilde{A}_j$ is normally
distributed with zero mean and standard deviation equal to 1.0,
then $A_j = \sigma_j \tilde{A}_j$ is normally distributed with
mean zero and standard deviation equal to $\sigma_j$. Similarly
for $B_j$. Eq.~\eqref{eq:8.1} may therefore be written in the
form,
\begin{equation}\label{eq:3.49}
{X}(t) = \sum_{j=1}^N  \sigma_j \, \big\{ \tilde{A}_j \, \cos
(\omega_j t) + \tilde{B}_j \, \sin (\omega_j t) \big\},
\end{equation}
where $\tilde{A}_j$ and $\tilde{B}_j$, $j = 1, \ldots, N$, is now
a set of independent, standard normally distributed variables.
There are computer programs that may be used to generate
independent outcomes of a standard, normally distributed variable.
It is seen that $2N$ outcomes are needed for this example.
Specifically, assume that $\tilde{a}_j$ and $\tilde{b}_j$, $j = 1,
\ldots, N$ are the obtained outcomes from such a program. The
corresponding realization $x(t)$ is then,
\begin{equation}\label{eq:3.50}
x(t) = \sum_{j=1}^N  \sigma_j \, \big\{ \tilde{a}_j \, \cos
(\omega_j t) + \tilde{b}_j \, \sin (\omega_j t) \big\}.
\end{equation}
This procedure can be repeated as many times as needed to
produce the requested sample size of realizations. It is, of
course, understood here that each realization is generated
independently of all others. This procedure is illustrated with an application to ocean waves in
Example~\ref{ex:51}.

The following representation is also
frequently used, especially when ergodic properties are desirable,

\begin{equation}\label{eq:8.6}
X(t) = \sum_{j=1}^N \, a_j \, \cos (\omega_j t + \Phi_j), 
\end{equation}
where $a_j$ and $\omega_j$, $j = 1, \ldots, N$ are positive
constants, and $\{\Phi_j\}_{j=1}^N$ is a set of independent random
variables that are uniformly distributed over $(0, 2 \pi)$. It may be shown that $X(t)$
is ergodic with respect to the mean value and the autocorrelation.


Generating realizations of this process is now quite
straight-forward. It is seen that, in fact, only outcomes of the
uniformly distributed phase angles $\Phi_j$, $j = 1, \ldots, N$
are needed. If $R_j$, $j=1, \ldots, N$, denote independent random
variables uniformly distributed on $(0,1)$, usually referred to as
random numbers, then clearly one may put $\Phi_j = 2 \pi R_j$.
Hence, by invoking a random number generator, usually available on
any computer, outcomes of $\Phi_j$, $j = 1, \ldots, N$ can be
easily generated. Let $r_j$, $j = 1, \ldots, N$ denote a set of
independent random number outcomes. A realization $x(t)$ of $X(t)$
is then,
\begin{equation}\label{eq:3.51}
x(t) = \sum_{j=1}^N \, a_j \, \cos (\omega_j t + 2 \pi r_j).
\end{equation}




\subsection{Variance spectra}\label{ssec:variancespectrum} \index{spectrum!variance} \index{variance spectrum}


Fourier analysis is routinely used to decompose the time
histories of load and response as sums or integrals of
$\cos(\cdot)$ and $\sin(\cdot)$ terms over the frequency domain.
It is known how the transfer function gives a direct
connection between the amplitudes at each frequency for the load
and response, provided the system is linear and time invariant. A
time history that is periodic can be decomposed as a sum over a
finite or countably infinite number of frequencies. If it is not
periodic, the decomposition must be expressed as an integral. It can be shown that this is possible only if the
time history dies out with time. Regarding the realizations of
stationary stochastic processes, they will not generally be
periodic. Nor will they decrease with time because the variance
is constant, and this is a measure of the fluctuations around the
mean value. A direct frequency decomposition of the realizations
of a stationary process is therefore not directly feasible. This
difficulty is circumvented by using the autocovariance
function of a stationary process. This function will generally approach zero when its time argument
increases.



Let $X(t)$ be a stationary process with autocovariance function
$C_X(\tau)$. Assume that $C_X(\tau) \to 0$ when $\tau \to \infty$
sufficiently rapidly so that $\int_{- \infty}^\infty |C_X(\tau)|
\, d \tau$ has a finite value, that is, $C_X(\tau)$ is assumed to be an integrable function.

The \emph{variance spectrum} $S_X(\omega)$ of $X(t)$ is defined as
the Fourier transform of $C_X(\tau)$ as follows,
\begin{equation}\label{eq:4.1}
S_X(\omega) = \frac{1}{2 \pi} \int_{- \infty}^\infty C_X(\tau) \,
e^{- \ii \omega \tau} \, d \tau.
\end{equation}

It was shown already in the 1930s that $C_X(\tau)$ and
$S_X(\omega)$ constituted a Fourier transform pair; that is,
$C_X(\tau)$ is given by the inverse Fourier transform as,
\begin{equation}\label{eq:4.2}
C_X(\tau) = \int_{- \infty}^\infty S_X(\omega) \, e^{\ii \omega
	\tau} \, d \omega.
\end{equation}

Equations~(\ref{eq:4.1}) and (\ref{eq:4.2}) are often called the
Wiener-Khintchine relations \index{Wiener-Khintchine relations} after the originators, cf. \cite{chat89}.

The quantity $S_X(\omega)$ is known by many names. A few names
that seem appropriate to mention here are \emph{energy spectrum}, \index{energy spectrum}
\emph{power spectral density}, \index{power spectral density} \emph{spectral density} \index{spectral density} or just 
\emph{spectrum}. \index{spectrum} The name \emph{variance spectrum}, ties in directly with the interpretation of $S_X(\omega)$. By
putting $\tau = 0$ in Eq.~(\ref{eq:4.2}), it is obtained that,
\begin{equation}\label{eq:4.3}
\sigma_X^2 = C_X(0) = \int_{- \infty}^\infty S_X(\omega) \, d
\omega.
\end{equation}
This equation shows that $S_X(\omega)$ can be interpreted as a
distribution of variance along the frequency axis, provided that
$S_X(\omega) \ge 0$, and in Section~\ref{sec:var-spec-real}, it is shown that this is
always the case. Negative frequencies are a mathematical
convenience and have no real physical content. If $S_X(\omega)$ is
to be interpreted as distribution of variance along the frequency
axis, one would therefore expect that $S_X(- \omega) =
S_X(\omega)$. This symmetry property of $S_X(\omega)$ can be shown
directly from Eq.~(\ref{eq:4.1}). It follows also from a rewriting of
the Wiener-Khintchine relations to real form, which is of interest
in itself. Because $C_X(\tau)$ is symmetric, it follows from
Eq.~(\ref{eq:4.1}) that,
\begin{align}\label{eq:4.4}
S_X(\omega) &= \frac{1}{2 \pi} \int_{0}^\infty C_X(\tau) \, e^{-
	\ii \omega \tau} \, d \tau + \frac{1}{2 \pi} \int_{0}^\infty C_X(-
\tau) \, e^{\ii \omega \tau} \, d \tau \nonumber \\
&= \frac{1}{\pi} \int_{0}^\infty C_X(\tau) \, \cos \omega \tau \,
d \tau \,,
\end{align}
where Euler's relation $e^{\ii x} = \cos x + \ii \sin x$ is used.
The last integral in the upper line is obtained by the change of
variable $\tau \to - \tau$, while the lower line follows from the
symmetry property $C_X(- \tau) = C_X(\tau)$.

From Eq.~(\ref{eq:4.4}), it follows that,
\begin{equation}\label{eq:4.5}
S_X(- \omega) = S_X( \omega),
\end{equation}
because $\cos(- \omega t) = \cos \omega t$. Eq.~(\ref{eq:4.2}) can
then be rewritten as,
\begin{equation}\label{eq:4.6}
C_X(\tau) = 2 \int_{0}^\infty S_X(\omega) \, \cos \omega \tau \, d
\omega.
\end{equation}

If $S_X( \omega)$ is a reasonably nice function, the integral on
the rhs of Eq.~(\ref{eq:4.6}) can be approximated by a finite sum,
viz.
\begin{equation}\label{eq:4.7}
C_X(\tau) = 2 \int_{0}^\infty S_X(\omega) \, \cos \omega \tau \, d
\omega \approx \sum_{j=1}^N 2 S_X(\omega_j) \, \Delta \omega \,
\cos \omega_j \tau,
\end{equation}
for a suitable choice of $\omega_1 < \ldots < \omega_N$, and
sufficiently small $\Delta \omega = (\omega_{N} -
\omega_1)/(N-1)$. Let us now recollect previous results and define a stationary process,
\begin{equation}\label{eq:4.8}
\tilde{X}(t) = \sum_{j=1}^N \, \big\{ A_j \, \cos \omega_j t + B_j
\, \sin \omega_j t \big\} \, ,
\end{equation}
where the random variables $A_j$ and $B_j$ satisfy the necessary
conditions for stationarity, cf. Eq.~\eqref{eq:8.1}. In addition, $\sigma_j^2 = \EE[A_j^2] = \EE[B_j^2] = 2
S_X(\omega_j) \Delta \omega$. Then, according to Eq.~\eqref{eq:4.3}, the autocovariance of $\tilde{X}(t)$ is given as,
\begin{equation}\label{eq:4.9}
C_{\tilde{X}}(\tau) = \sum_{j=1}^N 2 S_X(\omega_j) \, \Delta
\omega \, \cos \omega_j \tau.
\end{equation}
A stationary process $\tilde{X}(t)$ is thus constructed with approximately the same variance distribution, that is, variance spectrum, as ${X}(t)$. In a certain sense, $\tilde{X}(t)$
can be said to represent ${X}(t)$. What has just been described, is
a variant from a class of methods that is extensively used in
practice to generate realizations of a given stationary process.
To get a concrete realization of a process represented by
Eq.~(\ref{eq:4.8}), one has to generate outcomes of the random
variables that enter the sum.

On the basis of Eq.~(\ref{eq:4.5}) and the fact that negative
frequencies do not really have any physical meaning, it is common
practice in engineering to use the \emph{one-sided} variance
spectrum, which is denoted by $S_X^+(\omega)$ in this chapter, and defined as,
\begin{equation}\label{eq:4.10}
S_X^+(\omega) = \left\{ \begin{array}{l} 2 S_X(\omega) \, , \quad
\omega \ge 0, \\ 0\, , \quad \omega < 0.
\end{array} \right.
\end{equation}

The distribution of variance is thereby concentrated to positive
(physically realizable) frequencies. The variance expressed in
terms of the one-sided variance spectrum is clearly,
\begin{equation}\label{eq:4.11}
\sigma_X^2 = \int_{0}^\infty S_X^+(\omega) \, d \omega.
\end{equation}

A typical variance spectrum produced from measured data for the
wave elevation at a point on the sea surface at a location in the
North Sea is shown in Figure~\ref{fig:4.1}. The somewhat jagged
look is mostly due to low numerical resolution.


\begin{figure}[!hbt]
	\begin{center}
		%\setlength{\unitlength}{0.012500in}%
		\dofig{9cm}{Figures/fig4_1.pdf}
	\end{center}
	\caption{Typical "wave spectrum" from the North Sea.}
	\label{fig:4.1}
\end{figure}




\subsection{Units of variance spectra}\label{sec:uvs} \index{variance spectrum!units}

In Eq.~(\ref{eq:4.3}), the units of $S_X(\omega)$
are the square of the units of $X(t)$ divided by radians per
second. If $X(t)$ models the wave elevation at a location on the
ocean surface and is measured in meters, then the units of
$S_X(\omega)$ are m$^2$s/rad.

So far, the variance spectrum has only been considered as a function
of circular frequency $\omega$ with units rad/s. It is also quite
common to give the variance spectrum as a function of frequency
$f$ with units in Hertz (oscillations per second). The relation
$\omega = 2 \pi f$ implies a rescaling of the frequency axis when
$\omega$ is replaced by $f$. It is therefore not correct to believe
that $S_X(2 \pi f)$ would represent the variance spectrum as a
function of $f$; this would easily lead to the wrong variance. If $G_X^+(f)$ denotes the one-sided
variance spectrum as a function of $f$ in Hz, then to conserve the variance, it is needed to require that $G_X^+(f) df = S_X(\omega) d \omega$. 
This leads to the relation,

\begin{equation}\label{eq:4.12}
G_X^+(f ) = 2 \pi S_X^+(\omega).
\end{equation}

Exactly the same relation applies to two-sided spectra. A factor
of $2 \pi$ = 6.28 occurring erroneously in the variance can have a significant impact
on the results in some cases. Thus, care should be exercised when
adopting values of spectral moments by checking which kind of
frequency is used in the variance spectrum.




\subsection{Example - A realization of a wave process}\label{ex:51} \index{wave process!realization}


	
	In this example, the procedure for producing a
	realization of a stationary stochastic process described in
	Subsection~\ref{ssec:realizations} is illustrated. In particular, it will be shown
	how a realization of a stochastic process with a given variance
	spectrum can be generated. To be even more specific, it will be 
	assumed that the task is to produce an arbitrary realization of the
	wave elevation $X(t)$ at a given location in the North Sea with a
	recorded wave spectrum, as depicted in Figure~\ref{fig:4.1}. Thus,
	it is assumed that this wave elevation can be represented as a
	stationary stochastic process $X(t)$ with a one-sided variance
	spectrum $G_X^+(f)$, as shown in Figure~\ref{fig:4.1}. As mentioned
	in Section~\ref{sec:variancespectrum}, there are several
	alternative methods that can be used to generate realizations of a
	stochastic process. Our choice here is to use the method described
	in Section~\ref{sec:variancespectrum}, which amounts to
	approximating the process $X(t)$ with $\tilde{X}(t)$, and then
	generate realizations of $\tilde{X}(t)$ instead.
	
	To proceed, it is necessary to specify what kind of random
	variables $A_j$ and $B_j$ to use in Eq.~(\ref{eq:4.8}). The common
	choice in the case of wave processes on the open ocean is to
	assume that these variables are independent and normally
	distributed with zero mean value. To determine the standard
	deviation, it is necessary to decide on a suitable discretization of the
	frequency axis. For many practical purposes, it would be desirable
	to have a discretization that would give of the order of $10^3$
	frequencies in the frequency range where the waves have
	significant energy, e.g. from 0.27 to 0.99 Hz in this case, cf. the discussion at the end of the example.
	Because the objective here is to illustrate the procedure, 
 	a rather coarse frequency increment $\Delta f = 0.09$ Hz will do..
	That gives 8 frequencies denoted by $f_1, \ldots, f_8$ ($f_1 = 0.315$ Hz,
	$f_{j+1} = f_j + 0.09$ Hz) in the specified frequency range. The
	process $\tilde{X}(t)$ can then be written as,
	\begin{equation}\label{eq:4.13}
	\tilde{X}(t) = \sum_{j=1}^8 \, \big\{ A_j \, \cos (2 \pi f_j t) +
	B_j \, \sin (2 \pi f_j t) \big\}. 
	\end{equation}
	Figure~\ref{fig:4.4} shows the relevant part of the wave
	spectrum in Figure~\ref{fig:4.1} magnified along the frequency
	axis to clarify how the standard deviation of $A_j$ and $B_j$ is
	determined.
	
	\begin{figure}[!hbt]
		\begin{center}
			%\setlength{\unitlength}{0.012500in}%
			\dofig{7cm}{Figures/fig4_4.pdf}
		\end{center}
		\caption{Wave spectrum from the North Sea.} \label{fig:4.4}
	\end{figure}
	
	What remains to get an approximate realization of $X(t)$ is to
	generate outcomes of the random variables that enter
	$\tilde{X}(t)$. For this, an observation made in
	Subsection~\ref{ssec:realizations} is invoked: if $\tilde{A}_j$ is normally distributed with zero mean 
	and standard deviation equal to 1.0, then $A_j = \sigma_j
	\tilde{A}_j$ is normally distributed with zero mean and standard
	deviation equal to $\sigma_j$. Similarly for $B_j$. Because $\sigma_j =  \sqrt{G_X^+(f_j) \, \Delta f}$,  Eq.~\eqref{eq:4.13} may be written in the form,
	\begin{equation}\label{eq:4.21}
	\tilde{X}(t) = \sum_{j=1}^8  \sqrt{G_X^+(f_j) \, \Delta f} \,
	\big\{ \tilde{A}_j \, \cos (2 \pi f_j t) + \tilde{B}_j \, \sin (2
	\pi f_j t) \big\}, 
	\end{equation}
	where $\tilde{A}_j$ and $\tilde{B}_j$, $j = 1, \ldots, 8$, now
	constitute a set of independent, standard normally distributed variables.
	Computer programs are easily available for generating independent
	outcomes of a standard, normally distributed variable. It is seen
	that 16 is needed for the example.
	
	An alternative procedure is to use a table or computer program for
	generating (pseudo-)random numbers, which are uniformly
	distributed between 0 and 1. This can also be used by invoking the
	following result. If $\Phi(\cdot)$ denotes the distribution of a standard,
	normally distributed variable, and $R$ denotes a random variable
	that is uniformly distributed between 0 and 1, then the random
	variable $Z = \Phi^{-1}(R)$ is standard and normally distributed.
	By generating 16 independent outcomes of $R$: $r_1, \ldots,
	r_{16}$, then $z_1 = \Phi^{-1}(r_1), \ldots, z_{16} =
	\Phi^{-1}(r_{16})$ will be 16 independent outcomes of a standard,
	normally distributed variable. This procedure is used here, and
	the results are shown in Table~\ref{tab:4.1}, where $\tilde{a}_j =
	\Phi^{-1}(r_{2j-1})$ is an outcome of $\tilde{A}_j$ and $\tilde{b}_j = \Phi^{-1}(r_{2j})$ is an outcome of $\tilde{B}_j$.
	
	
	\begin{table}[!htb] 
		\caption{A table of the numerical information needed to produce a realization of the stochastic process given by Eq.~\eqref{eq:4.21}.}\label{tab:4.1}
		\begin{center}
			\begin{tabular}{|c|c|c|c|c|c|c|}
				\hline
				$j$ & $f_j$[Hz] & $\sqrt{G_X^+(f_j) \, \Delta f}$ [cm] & $r_{2j-1}$ & $r_{2j}$ & $\tilde{a}_j$ & $\tilde{b}_j$ \\
				\hline \hline
				$1$ & $0.315$ & $1.90$ & 0.10097 & 0.32533 & -1.2760 & - 0.4529 \\
				\hline
				$2$ & $0.405$ & $5.45$ & $0.76520$ & 0.13586 & 0.7232 & - 1.0991  \\
				\hline
				
				3 & 0.495 & 6.57 & 0.34673 & 0.54876 & - 0.3942 & 0.1225 \\
				\hline
				
				4 & 0.585 & 3.15 & 0.80959 & 0.09117 & 0.8764 & - 1.3335 \\
				\hline
				
				5 & 0.675 & 2.23 & 0.39292 & 0.74945 & - 0.2717 & 0.6728 \\
				\hline
				
				6 & 0.765 & 1.56 & 0.37542 & 0.04805 & - 0.3175 & - 1.6641 \\
				\hline
				
				7 & 0.855 & 1.44 & 0.64894 & 0.74296 & 0.3825 & 0.6525 \\
				\hline
				
				8 & 0.945 & 1.41 & 0.24805 & 0.24037 & - 0.6807 & - 0.7051 \\
				\hline
				
				\hline
			\end{tabular}
		\end{center}
	\end{table}
	
	
	A piece of the corresponding realization is shown in
	Figure~\ref{fig:4.5}, and one may get a similar impression as when
	observing irregular seas out on the oceans. In practice, there is
	often a need to generate many realizations to perform 
	statistical analyses. It is then necessary to repeat the procedure
	just described the requisite number of times, and for each realization, a new set of outcomes independent of the previous ones are chosen.
	
	
	\begin{figure}[!hbt]
		\begin{center}
			%\setlength{\unitlength}{0.012500in}%
			\dofig{6cm}{Figures/fig4_5.pdf}
		\end{center}
		\caption{A realization of the wave elevation.} \label{fig:4.5}
	\end{figure}
	
	
	If a piece of the realization that has been generated, was shown
	with a duration which was twice as long, it would have become clear that the wave
	pattern is repeating itself. This is due to the way it was
	constructed, which indeed makes it periodic. The period is
	determined by the greatest common divisor of the frequency
	increment and the initial frequency of the discretization of the
	frequency range that is chosen. In this case, the
	period becomes 1/0.045 = 22.2~s. The practical consequence
	of this is that one must choose a discretization that is in
	correspondence with the required length of a realization. It may
	be worth mentioning that for most practically relevant
	discretizations, one may say that a corresponding realization will
	have a period that can be assumed to be given as $1/ \Delta f$.
	One does not fully avoid the problem related to periodicity by choosing
	an (almost) irrational ratio between the starting frequency and
	the frequency increment, or other ''smart'' tricks such as choosing the 
	frequencies randomly within each subinterval of the discretization.
	
	Another point worth noting is that the process $\tilde{X}(t)$ is
	not ergodic. If it is desirable to ensure this property, one may
	use the method described in Subsection~\ref{ssec:realizations}. However, the
	difference in the practical results obtained by using this method
	versus the one described here, is usually rather small if
	the discretization is properly done.
	


\subsection{The variance spectrum directly from the realizations}\label{sec:var-spec-real} \index{variance spectrum!from realizations}


When the variance spectrum was defined, it was mentioned that a
realization $x(t)$ of a stationary process $X(t)$ does not have a
Fourier transform because it does not decrease toward zero for
large $t$. Intuitively, one would nevertheless expect that most of
the information regarding the frequency content of $x(t)$ should
be contained in a finite section of $x(t)$ if this section is
large enough. A section of $x(t)$ can be defined as follows,
\begin{equation}\label{eq:4.33}
x_T(t) = \left\{ \begin{array}{l} x(t) \, , \quad 0 \le t \le T,
\\
0\, , \quad \mbox{elsewhere}.
\end{array} \right.
\end{equation}
Because $x_T(t)$ is zero outside a finite interval, it has a Fourier
transform, viz.
\begin{equation}\label{eq:4.34}
X_T(\omega) = \frac{1}{2 \pi} \int_{- \infty}^\infty  x_T(t) \,
e^{ - \ii \omega t} \, d t = \frac{1}{2 \pi} \int_{0}^T  x(t) \,
e^{ - \ii \omega t} \, d t.
\end{equation}
It would seem natural to expect that there is a connection between
$X_T(\omega)$ and the variance spectrum $S_X(\omega)$. And it
turns out that the connection is, in fact, quite simple. It is written in the following way: assume that $x_j(t)$, $j = 1, 2, \ldots$ are realizations of a stationary process. Then the
following equation applies,
\begin{equation}\label{eq:4.35}
S_X(\omega) = \lim_{T \to \infty} \lim_{N \to \infty} \frac{2
	\pi}{T N} \sum_{j=1}^N | X_{j,T}(\omega)|^2,
\end{equation}
where $ X_{j,T}(\omega)$ denotes the Fourier transform of
$x_{j,T}(t)$, which equals $x_j(t)$ for $0 \le t \le T$ and zero
elsewhere. From Eq.~(\ref{eq:4.35}), it is also immediately seen that
$S_X(\omega) \ge 0$.

Equation~(\ref{eq:4.35}) is based on the availability of an ensemble of
realizations. As previously discussed, there are many situations
where only one realization is available. By assuming that the
process is ergodic, it can be shown that the variance spectrum can
be determined in the following way. Let $x(t)$ denote a
realization of the assumed ergodic process $X(t)$. Define a set of
truncated Fourier transforms $X_{(j), T}(\omega)$ over the intervals $\big((j-1)T, jT\big)$, $j = 1, 2, \ldots$, as follows,
\begin{equation}\label{eq:4.36}
X_{(j),T}(\omega) = \frac{1}{2 \pi} \int_{(j-1)T}^{jT} x(t) \, e^{
	- \ii \omega t} \, d t.
\end{equation}
It is then obtained that,
\begin{equation}\label{eq:4.37}
S_X(\omega) = \lim_{T \to \infty} \lim_{N \to \infty} \frac{2
	\pi}{T N} \sum_{j=1}^N | X_{(j),T}(\omega)|^2.
\end{equation}
Equation~(\ref{eq:4.35}) or (\ref{eq:4.37}) can be regarded as the
basis for the \emph{fast Fourier transform} (FFT) algorithm for
calculating spectra. This method, developed around 1965, has
assumed a dominating position among numerical methods for
calculating Fourier transforms. This is primarily due to the fact
that the method is much faster than traditional methods. An
extensive discussion of the FFT method for calculating Fourier
transforms is given by \cite{newl91}.





\section{Monte Carlo Simulation of Load and Response} \label{sec:sstr} \index{Monte Carlo simulation}

The linear or linearized equations of motion for marine structures considered in
this book and their solutions are discussed in \cite{namo13}. For such equations, the comments of the previous section apply. 
An example of a nonlinear dynamic model that is often adopted for an offshore structure, can be written in the following form,
\begin{equation}\label{eq:12.1}
\bm{M}\, \ddot{\bm{X}}(t) + \bm{C}\, \dot{\bm{X}}(t) +
\bm{f}\big(\bm{X}(t), \dot{\bm{X}}(t), t \big) = \bm{P}(t)\,,
\end{equation}
where $\bm{M}$ and $\bm{C}$ are suitable mass and damping matrices, respectively; $\bm{f}(\cdot, \cdot) = \big(
f^{kl}(\cdot, \cdot)\big)$, $k,l = 1, \ldots, n$, is a nonlinear
matrix function; $\bm{P}(t) = (P_1(t), \ldots, P_n(t))'$ denotes a
stochastic loading process; and $\bm{X}(t) = (X_1(t),$ $ \ldots, X_n(t))'$
is the corresponding response process.

The Monte Carlo method applied to such a system would consist of
generating a statistical sample of specified size $N$, say, of
response time histories by first generating a sample of time
histories of the same size of the loading process, and then solve
Eq.~(\ref{eq:12.1}) for each of the load time histories in the
sample using methods of numerical integration of dynamical system equations, cf. \cite{arml91,namo13}. When the desired
sample of response time histories is produced, the
statistical analysis of the response may then proceed as discussed
in the remaining part of this chapter.



\section{Sample Statistics of Simulated Response} \label{sec:sssr} \index{sample statistics}

An important element for graphic representation of sampled data
from a statistical population for which the underlying
distribution function is unknown is the so-called plotting
position formula.\index{plotting position formula} It is based on the notion of order statistics.\index{order statistics}

Let us start by assuming that $X$ is a continuous random variable
with a distribution $F(x)$ and a probability density $f(x)$. The given sample of
independent observations $x_1, x_2, \ldots, x_n$ is now ordered in
an increasing sequence $x_{(1)} \le x_{(2)} \le \ldots \le
x_{(n)}$. The random variable $X_{(m)}$ corresponding to $x_{(m)}$
is called the $m$th order statistic, $m = 1, \ldots, n$. The probability density
$f_m(x)$ of $X_{(m)}$ follows from the observation that
$X_{(m)}=x$ implies the event that there are $m-1$ outcomes of $X$
with values less than (or equal to) $x$, and $n-m$ outcomes with
values greater than $x$. According to the binomial distribution,
the probability of this event equals $n !/\big((m-1)! (n-m)!\big)
\, F(x)^{m-1} \big(1 - F(x)\big)^{n-m}$. Hence, it is obtained that \citep{cabe02},
\begin{equation}\label{eq:11.4.1}
f_m(x) = \frac{n !}{(m-1)! (n-m)!}\, F(x)^{m-1} \big(1 - F(x)
\big)^{n-m} \, f(x).
\end{equation}

It is now required to calculate $\EE[ F(X_{(m)})]$. This is given by,
\begin{align}\label{eq:11.4.2}
 &\EE[ F(X_{(m)}) ] = \int_x F(x) f_m(x)\, dx =
m {n \choose m}\int_x F(x)^m \big(1 - F(x) \big)^{n-m} \, f(x)
\,dx \nonumber\\
&= m {n \choose m} \int_0^1 F^m (1 - F)^{n-m} \, dF = m {n \choose
m} \frac{m!(n-m)!}{(n+1)!} = \frac{m}{n+1}.
\end{align}


Similarly, the variance of $F(X_{(m)})$ is calculated to be,
\begin{align}\label{eq:11.4.3}
\Var[F(X_{(m)}] &= \EE[ F(X_{(m)})^2 ] - \EE[ F(X_{(m)})]^2 \nonumber\\
&= \frac{m (m+1)}{(n+1)(n+2)} - \frac{m^2}{(n+1)^2} = \frac{m (n+1
- m)}{(n+1)^2 (n+2)}.
\end{align}


The results of Eqs.~(\ref{eq:11.4.2}) and (\ref{eq:11.4.3}) are
useful because they provide a means of plotting the sample of
observations $x_1, x_2, \ldots, x_n$ in order to estimate the distribution
$F(x)$ empirically. Equation~(\ref{eq:11.4.2}) states that the expected
value of the distribution function evaluated at the observation of
order $m$ is equal to $m/(n+1)$. This result suggests that an
optimal plotting strategy is obtained by plotting the points
$\big(x_{(m)}, m/(n+1)\big)$, $m=1, \ldots, n$.
Equation~(\ref{eq:11.4.3}) provides information on the variance of
$F(X_{(m)})$, that is, the ordinate of the plotting point. Due to
the symmetry of the expression, it attains its
maximum $1/[4(n+2)]$ at the median and decreases symmetrically to
$n/[(n+1)^2(n+2)]$ toward the ends of the interval $[0, 1]$.
Because the distribution function itself is not known, these results
are called distribution-free\index{distribution-free} results.

The preceding results are complemented by calculating the
probability that the $m$th observation is not exceeded by a future
observation. Now, the conditional probability that a single
observation will not exceed $X_{(m)}$ given that $X_{(m)} = x$ is
equal to $F(x)$. The corresponding unconditional probability $p_m$, say, is then obtained by using the law of total probability,
which gives exactly the same result as given by
Eq.~(\ref{eq:11.4.2}), that is,
\begin{align}\label{eq:11.4.4}
p_m = \int_x F(x) f_m(x)\, dx = \frac{m}{n+1}.
\end{align}
This result shows that a new observation of the continuous random
variable $X$ has equal probability of assuming a value in any of
the $n+1$ intervals defined by the previous $n$ observations. This
lends further support to the optimality of the plotting position
formula expressed by Eq.~\eqref{eq:11.4.2}.


It may be pointed out that several alternative plotting position
formulas for specific classes of distributions have been suggested
over the years, which primarily aim at correcting for sample bias.
A discussion of this topic is not pursued here, but the reader is rather referred to the literature. A discussion of particular interest related to the estimation of return periods can be found in \cite{makk06,makk08}, where it is argued for the use of Eq.~\eqref{eq:11.4.2}. 

A useful diagnostic tool to check the accuracy of an assumed
statistical distribution $F$ for the observed data is obtained
by comparing the fitted distribution $\hat{F}$ with the data on
a quantile plot\index{quantile plot} (QQ-plot)\index{QQ-plot} or a probability plot\index{probability plot} (PP-plot). \index{PP-plot}
Assuming that $\hat{F}$ is strictly increasing and continuous, the
QQ-plot is obtained by comparing the ordered data with the
corresponding quantiles of the fitted distribution by plotting
\begin{align}\label{eq:11.4.5}
\Big(\hat{F}^{-1}\big(\frac{m}{n+1} \big), x_{(m)} \Big) \,.
\end{align}
The name \emph{QQ-plot} derives from the fact that both
$\hat{F}^{-1}\big(\frac{m}{n+1} \big)$ and $x_{(m)}$ are estimates
of the $m/(n+1)$th quantile of $F$. If $F$ is a good choice for
the distribution of the data, the QQ-plot should be close to the
straight line of slope 1 passing through the origin.
Alternatively, the fit of $\hat{F}$ to the data can be checked by
the PP-plot, which is obtained by plotting the points
\begin{align}\label{eq:11.4.6}
\Big(\hat{F}\big(x_{(m)} \big), \frac{m}{n+1} \Big) \,.
\end{align}
A good fit is again demonstrated if the plotted graph is
approximately a straight line. The main difference between the two
plots is that the QQ-plot gives a more clear impression of the fit
of the tail data, which may be of particular significance for
extreme value statistics.

It is worth noting that for a range of distribution functions
QQ-plots can be constructed without having to estimate
distribution parameters. This typically applies to distributions
characterized by a scale and a location parameter, e.g. the normal
distribution. In such cases the intercept of the line fitted to
the QQ-plot would represent location, while the slope represents
scale. For example, a QQ-plot for a normal distribution can be
achieved by plotting $\big(\Phi^{-1}(m/(n+1)), x_{(m)} \big)$.



\section{Latin Hypercube Sampling} \label{sec:lhcs} \index{Latin hypercube sampling}

Latin hypercube sampling (LHS)\index{LHS} \citep{mccobe79} is a method for
effectively reducing the sample size for Monte Carlo simulations
of stochastic response processes that depend on many random
parameters. In cases where the parameters of a dynamic model, such as
mass, damping, and stiffness are modeled as random variables, the
response process will then also depend on these random variables.
Specifically, let us assume that the model depends on the random
parameters $Y_1, \ldots, Y_m$. To highlight the dependence on
these parameters, a random response process $X(t)$ of this model
may then be written as $X(t) = X(t;Y_1, \ldots, Y_m)$. Because the
external loading is often modeled as a stochastic process, the
response $X(t;y_1, \ldots, y_m)$ becomes a stochastic process for
each sample $y_1, \ldots, y_m$ of the random parameters. If $m$ is
not small, then the number of samples needed to provide good
sample statistics for the response process may become huge
if no consideration is made on how to effectively represent the
statistical variability of the random parameters. LHS is a very
good and simple method for this purpose.

LHS starts by selecting $k$ different values from each of the $m$
random variables in the following manner. The interval $(0,1)$ is
divided into $k$ equally long intervals $I_j = ((j-1)/k, j/k)$, $j
= 1, \ldots, k$. Let $F_i(y)$ denote the distribution function of
$Y_i$, $i = 1, \ldots, m$. For each $i$, $k$ independent outcomes
$u_1, \ldots, u_k$ of the random number $U$, which is uniformly
distributed on (0,1), are produced. The resulting ordered sample
for $Y_i$ is $y_{i,1} < \ldots < y_{i,k}$ where $y_{i,j} =
F_i^{-1}((j-1)/k + u_j/k)$. Note that $(j-1)/k + u_j/k$ is nothing
but an outcome of a random number uniformly distributed on the
interval $((j-1)/k, j/k)$. A Latin hypercube (LH) presample for $Y_i$
is then obtained as $y_{i,r_1}, \ldots, y_{i,r_k}$, where $r_1,
\ldots, r_k$ is a random reordering of $1, \ldots, k$. Finally, an
LH sample for $Y_1, \ldots, Y_m$, which will also be
of size $k$, is now represented by the $m \times k$ array or
matrix $(y_{i,r_j})$, where each column is an element in the LH sample.

It is tacitly assumed that the random parameters are
independent. If this is not the case, LHS can also deal with
correlated parameters. Standard statistical software packages
usually offer LHS as an optional sampling technique.



\section{Estimation of Extreme Response}\label{sec:eer} \index{extreme response!estimation}

The view that estimation of extreme values by Monte Carlo methods is
generally prohibitive in terms of computer time, is a truth in
contention. In fact, with present-day computational power, it is
possible to perform simulations on a scale that allows
estimation of extreme structural response for a range of problems.
Admittedly, it is not difficult to describe a dynamic model for
which the required simulation time for direct estimation of
extreme response would be beyond any acceptable bounds, but the
systems for which this is the case are steadily diminishing in
step with the development of increasingly powerful computers. This
being the case, it is good reason to believe that Monte Carlo
simulation-based methods will become much more frequently used
even for extreme value estimation than what is the case today. It
is therefore considered appropriate to discuss such methods here.
Note that the methods described in this section are equally
applicable to measured response time histories, obtained, for
instance, from real life experiments. In this section, only the Gumbel and the point process methods are applied. For simplicity, only the case of positive extreme values will be discussed. The necessary
modifications to deal with other situations are usually obvious.






\subsection{The Gumbel method}\label{ssec:gm}\index{Gumbel method}\index{extreme value estimation!Gumbel method}

It was pointed out in Chapter~\ref{chap:acer} that for response
processes relevant for many engineering structures, the appropriate extreme
value distribution would almost always be the Gumbel distribution.
Therefore, let us assume that this is indeed the case for the
response process $X(t)$, which can be simulated by a suitable
procedure. Now, assume that $N$ independent response time
histories, each of duration $T$, have been simulated for a given
environmental condition. For the Gumbel method, the extreme
response is then identified for each time series or block of data. These extreme
value data are assumed to be Gumbel distributed, and plotting the
obtained data set of extreme values using a Gumbel probability plot
should then ideally result in a straight line. In practice, one
cannot expect this to happen, but on the premise that the data
follow a Gumbel distribution, a straight line can be fitted to the
data. Due to its simplicity, a popular method for fitting this straight line is the
method of moments, cf. Section~\ref{sec:pegev}. That is, writing the Gumbel distribution of the
extreme value $M(T)$ as,
\begin{align}\label{eq:11.7.1}
\Prob( M(T) \le \xi ) = \exp \left\{ - \exp \big( - a ( \xi - b)
\big) \right\},
\end{align}
it was shown that the parameters $a > 0$ and $b$ are related to the
mean value $m_M$ and standard deviation $\sigma_M$ of $M(T)$ as
follows: $b = m_M - 0.5772\, a^{-1}$ and $a = 1.2826 / \sigma_M$
\citep{bury75}. The estimates of $m_M$ and $\sigma_M$ obtained
from the available sample therefore provide estimates of $a$ and
$b$, which leads to the fitted Gumbel distribution by the method of moments.

Typically, a specified fractile value of the fitted Gumbel
distribution is then extracted and used in a design consideration.
To be specific, let us assume that the requested fractile value is
the $100(1-\alpha)$\% fractile, where $\alpha$ is usually a small
number, for example $\alpha = 0.1$. To quantify the uncertainty
associated with the obtained $100(1-\alpha)$\% fractile value
based on a sample of size $N$, the 95\% confidence interval of
this value is often used. A good estimate of this confidence
interval can be obtained by using a parametric bootstrapping
method \citep{efti93,dahi97}, cf. Section~\ref{sec:bootstrapping!GEV}. 
In our context, this simply means that the initial sample of $N$ extreme values is assumed to have
been generated from an underlying Gumbel distribution, whose
parameters are, of course, unknown. If this Gumbel distribution
had been known, it could have been used to generate a large number
of (independent) samples of size $N$. For each sample, a new
Gumbel distribution would be fitted and the corresponding
$100(1-\alpha)$\% fractile value identified. If the number of
samples had been large enough, an accurate estimate of
the 95\% confidence interval on the $100(1-\alpha)$\% fractile
value based on a sample of size $N$ could be found. Because the true
parameter values of the underlying Gumbel distribution are
unknown, they are replaced by the estimated values obtained from
the initial sample. This fitted Gumbel distribution is then used
as previously described to provide an approximate 95\% confidence
interval. Note that the assumption that the initial $N$ extreme
values are actually generated with good approximation from a
Gumbel distribution cannot be easily verified, which is
a drawback of this method. As has been pointed out, compared with the POT method, the
Gumbel method would also seem to use much less of the information
available in the data. This may explain why the POT method has
become increasingly popular over the past years, but the Gumbel
method is still widely used in practice.




\subsection{The point process method}\label{ssec:ngm}\index{point process method}\index{extreme value estimation!point process method}

It is known from Eq.~(\ref{eq:6.30}) that a reasonably good approximation of the
distribution of the extreme value $M(T)$ is obtained from the formula,
\begin{align}\label{eq:11.3.1}
\Prob( M(T) \le \xi ) = \exp \{ - \nu^+(\xi) \, T \},
\end{align}
where $\nu^+(\xi)$ denotes the mean upcrossing rate of a stationary process $X(t)$.

The method to be discussed here, relies on this particular
approximation. This implies that the mean upcrossing rate needs to
be estimated from the simulated time series. Assuming the
so-called ergodic mean value property, it is obtained that
\begin{align}\label{eq:11.3.2}
\nu^+(\xi) = \lim_{T \to \infty} \frac{1}{T} \, n^+(\xi;T)\,,
\end{align}
where $n^+(\xi;T)$ denotes a realization of $N^+(\xi;T)$; that is,
$n^+(\xi;T)$ denotes the counted number of upcrossings during time
$T$ from a particular simulated time history. In practice, $k$
time histories of a suitable length $T_0$, say, are provided by
simulation. The appropriate ergodic mean value estimate of
$\nu^+(\xi)$ is then,
\begin{align}\label{eq:11.3.3}
\hat{\nu}^+(\xi) = \frac{1}{k \,T_0} \sum_{j=1}^k  \,
n_j^+(\xi;T_0)\,,
\end{align}
where $n_j^+(\xi;T_0)$ denotes the counted number of upcrossings
of the level $\xi$ from time history no. $j$. This will often be the
chosen approach to the estimation of the mean upcrossing rate.

For a suitable number $k$, e.g. $k \ge 20-30$, and provided that
$T_0$ is sufficiently large, a fair approximation of the 95\%
confidence interval (CI$_{0.95}$) for the value $\nu^+(\xi)$ can
be obtained as,
\begin{equation}\label{eq:11.3.4}
\mbox{CI}_{0.95}(\xi)= \Big(\hat{\nu}^+(\xi) - 1.96 \,
\frac{\hat{s}(\xi)}{\sqrt{k}} \,, \: \hat{\nu}^+(\xi) + 1.96 \,
\frac{\hat{s}(\xi)}{\sqrt{k}} \,\Big)\,,
\end{equation}
where the empirical standard deviation $\hat{s}(\xi)$ is given as,
\begin{align}\label{eq:11.3.5}
\hat{s}(\xi)^2 = \frac{1}{k-1} \sum_{j=1}^k  \, \Big(
         \frac{n_j^+(\xi;T_0)}{T_0} - \hat{\nu}^+(\xi) \, \Big)^2.
\end{align}
Note that $k$ and $T_0$ may not necessarily be the number and
length of the actually simulated response time series. Rather,
they were chosen to optimize the estimate of
Eq.~(\ref{eq:11.3.5}). If, initially, $\tilde{k}$ time series of
length $\tilde{T}$ are simulated, then $k = \tilde{k} \, k_0$ and
$\tilde{T} = k_0 \, T_0$. That is, each initial time series of
length $\tilde{T}$ was divided into $k_0$ time series of
length $T_0$, assuming, of course, that $\tilde{T}$ is large
enough to allow for this in an acceptable way. The consistency of
the estimates obtained by Eq.~(\ref{eq:11.3.5}) can be checked by
the observation that $\Var[N^+(\xi;t)] = \nu^+(\xi)\,t$ when
$N^+(\xi;t)$ becomes a Poisson random variable, which by
assumption occurs for large values of $\xi$. This leads to the
equation,
\begin{align}\label{eq:11.3.6}
{s}(\xi)^2 = \frac{1}{k} \, \Var\left[  \sum_{j=1}^k \,
\frac{N_j^+(\xi;T_0)}{T_0} \right] = \frac{\nu^+(\xi)}{T_0},
\end{align}
where $\{N_1^+(\xi;T_0), \ldots, N_k^+(\xi;T_0) \}$ denotes a
random sample with a possible outcome $\{n_1^+(\xi;T_0), \ldots,
n_k^+(\xi;T_0) \}$. Hence, $\hat{s}(\xi)^2/k \approx
\hat{\nu}^+(\xi)/ (k \, T_0)$. Because this last relation is
consistent with the adopted assumptions (for large $\xi$), it
could have been used as the empirical estimate of the variance in
the first place. It is also insensitive to the blocking of data
discussed previously because $k \, T_0 = \tilde{k}\tilde{T}$. However,
the accuracy of this approach may be poor for small to moderate
values of $\xi$, where the Poisson assumption about the upcrossing
events may fail. In contrast, the advantage of
Eq.~(\ref{eq:11.3.5}) is that it does not rely on any specific
assumptions about the statistical distributions involved.

The idea underlying the development of the approach described here is based on the observation that
for dynamic models relevant for most engineering structures, the mean
$\xi$-upcrossing rate as a function of the level $\xi$ is 
highly regular in a particular way \citep{naga08jem}. As is shown
later, the mean upcrossing rate tail, say, for $\xi \ge \xi_0$,
behaves similarly to $\exp\{- a( \xi - b)^c \}$ ($\xi \ge
\xi_0$), where $a > 0$, $b \le \xi_0$, and $c > 0$ are suitable constants.
Hence, as discussed in detail by \cite{naga08jem}, it may be
assumed that,
\begin{align}\label{eq:11.3.7}
\nu^+(\xi) \approx q(\xi) \, \exp\{- a(\xi - b)^c \} \, , \: \xi
\ge \xi_0,
\end{align}
where the function $q(\xi)$ is slowly varying compared with the
exponential function $\exp\{-  a(\xi - b)^c \}$.
Equation~(\ref{eq:11.3.7}) can be rewritten as,
\begin{align}\label{eq:11.3.8}
\ln \big| \ln \big(\nu^+(\xi)/q(\xi) \big)\big| \approx c\,
\ln(\xi - b) + \ln a \, , \: \xi \ge \xi_0 \,.
\end{align}

It follows that by plotting $\ln \big| \ln
\big(\nu^+(\xi)/q(\xi) \big)\big|$ versus $\ln (\xi - b)$, it is
expected that an almost perfectly linear tail behavior will be
obtained. Now, as it turns out, the function $q(\xi)$ can be
largely considered as a constant $q$, say, for tail values of
$\xi$. This suggests using a method for identifying the parameters
$q$ and $b$ by optimizing the linear fit in the tail. When this
is achieved, the corresponding values of $a$ and $c$ can
then be extracted from the plot. This is discussed at some length
in \cite{naga08jem}. A plot of $\ln \big| \ln \big(\nu^+(\xi)/q
\big)\big|$ versus $\ln(\xi - b)$ for optimal parameters $b$ and
$q$ will be referred to as an optimal transformed plot. Examples of such plots will be given for some of the examples to follow, mostly for the purpose of demonstrating the validity of assumptions. An alternative, more extensive method for optimizing the fit to the data is described in Chapter~\ref{chap:acer}.  




In engineering applications, it is quite common to assume that the
observed extreme value response data do follow a Gumbel
distribution, cf. Subsection~\ref{ssec:gm}. The problem with this
approach is that classical extreme value theory cannot be used to
decide to what extent the asymptotic distribution is actually
valid for a given set of extreme value data. Note that the
asymptotic Gumbel distribution given by Eq.~(\ref{eq:11.7.1})
corresponds to an asymptotic upcrossing rate that is purely
exponential, that is, with $c = 1$ in Eq.~(\ref{eq:11.3.7}).
Hence, by adopting a much more general class of functions, with
the purely exponential functions as a subclass, to represent the
upcrossing rate, as done here, the ability to capture
subasymptotic behavior is greatly enhanced. By this, the
necessity to adopt a strictly asymptotic extreme value
distribution of questionable validity is avoided. However, note that for any $c > 0$, the corresponding extreme value model will be asymptotically Gumbel.

Note also that the so-called Weibull method for extreme value prediction, which is based on the assumption that the local peak values follow a three parameter Weibull distribution is, in fact, basically a prejudiced version of the point process  method in the sense that the parameter $q$ is a priori given the value 1. 

In cases where the approximation implied by Eq.~\eqref{eq:11.3.1} may be questioned, the ACER method should be applied.






\subsection{A comparison of methods}\label{ssec:cm}

In this subsection, the performance of the point process
method\index{point process method} is compared with that of the Gumbel method.\index{Gumbel method} This is done for the particular case of the horizontal deck response of a jacket
structure installed on the Kvitebj{\o}rn field in the North Sea. For this kind of response process, the point process approach is very accurate because a plot of the ACER functions shows that beyond the second ACER function, which corresponds to the upcrossing rate for suitably sampled time series, there are no dependence effects that need to be accounted for. 

Figure~\ref{fig:11.7.f1} depicts the Kvitebj{\o}rn jacket platform\index{Kvitebj{\o}rn jacket platform}
with the superstructure removed together with the corresponding three-dimensional
computer model used for the Monte Carlo simulations, see \cite{nagaha07oe} for details.



\begin{figure}[!tbh]
\begin{center}
\begin{tabular}{ccc}
\dofig{2.5cm}{Figures/Kvitebjorn_ss.pdf} & {\hspace{1.5cm}} &
\dofig{2.5cm}{Figures/Kvitebjorn_cm.pdf} \\
\end{tabular}

\caption{Left: Sketch of the Kvitebj{\o}rn platform with
the superstructure removed. Right: Computer model of
the Kvitebj{\o}rn platform \citep{kahabasp01}.}\label{fig:11.7.f1}
\end{center}
\end{figure}




\begin{table}[!hbt]
\begin{center}
\caption{Representative sea states. } \label{tab:11.7.1}
\begin{tabular}{c c}
\hline
$H_s$ (m) & $T_p$ (s)\\
\hline
12.0 & 12.0\\%3
14.7 & 16.5\\%1 218
\hline
\end{tabular}
\end{center}
\end{table}


For the simulations discussed here, two long-crested sea states
described by a JONSWAP wave spectrum as listed in
Table~\ref{tab:11.7.1}, were used. Twenty independent response time
histories, each of 3 hours' duration, were simulated for each sea
state. For the Gumbel method, the extreme horizontal deck response
in the wave direction is identified for each time series. These
extreme value data are then assumed to be Gumbel distributed, and
plotting each data set as a Gumbel probability plot
results in Figures~\ref{fig:11.7.f2b} and \ref{fig:11.7.f2c}.
Specifically, the observed 3-hour extremes $M_k$ are plotted versus
$-\ln \big(\ln(21/k)\big)$, for $k=1, \ldots, 20$. The fitted
straight line in each figure, which represents the fitted Gumbel
distribution, is based on the moment estimation method, cf.
Subsection~\ref{ssec:gm}.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Gumbel20_H12T12.pdf}
\end{center}
\caption{Empirical Gumbel plot of the 20 simulated 3-hour extremes of the horizontal deck displacement for
the sea state $H_s$ = 12~m and $T_p$ = 12~s together with the
fitted Gumbel distribution ($---$)} \label{fig:11.7.f2b}
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Gumbel20_H14_7T16_5.pdf}
\end{center}
\caption{Empirical Gumbel plot of the 20 simulated 3-hour extremes of the horizontal deck displacement for
the sea state $H_s$ = 14.7~m and $T_p$ = 16.5~s together with the
fitted Gumbel distribution ($---$)} \label{fig:11.7.f2c}
\end{figure}


The 90\% fractile value $L_G$ of the fitted Gumbel distribution is
identified and shown in each figure. Table~\ref{tab:11.7.2} lists the obtained 90\% fractile values. To quantify the
uncertainty associated with the obtained 90\% fractile value based
on a sample of size 20, the 95\% confidence interval (CI$_{0.95}$)
of this value is used. A good estimate of this confidence
interval can be obtained by using a parametric bootstrapping
method \citep{efti93,dahi97}. In our context, this simply means
that the initial sample of 20 extreme values is assumed to have
been generated from an underlying Gumbel distribution, whose
parameters are, of course, unknown. If this Gumbel distribution
had been known, it could have been used to generate many samples
of size 20. For each sample, a new Gumbel distribution would be
fitted, and the corresponding 90\% fractile value identified. If
the number of samples had been large enough, an accurate
estimate of the 95\% confidence interval on the 90\% fractile
value based on a sample of size 20 could be found. Because the true
parameter values of the underlying Gumbel distribution are
unknown, they are replaced by the estimated values obtained from
the initial sample. This fitted Gumbel distribution is then used
as previously described to provide an approximate 95\% confidence
interval. Note that the assumption that the initial 20 extreme
values are actually generated from a Gumbel distribution is quite
accurate in this case, as discussed later.


Invoking the parametric bootstrap,\index{bootstrapping!parametric} the 95\% confidence interval is
estimated for each case based on 100,000 bootstrap samples. The
obtained results are listed in Table~\ref{tab:11.7.2}. This way of
estimating the 90\% fractile value of the 3-hour extreme value
distribution is referred to as the Gumbel method. The
empirical densities obtained for the predicted 90\% fractile values
with the CI$_{0.95}$ indicated are shown in
Figures~\ref{fig:11.7.f2e} and \ref{fig:11.7.f2f}.


\begin{table}[!hbt]
\begin{center}
\caption{90\% fractile values $L_G$ of the fitted Gumbel distributions }\label{tab:11.7.2}
\begin{tabular}{c c c c}
\hline
$H_s$ (m) & $T_p$ (s) &  $L_G$ (m) & CI$_{0.95}$ \\
\hline
12.0 & 12.0 & 0.47 & (0.40, 0.54)\\
14.7 & 16.5 &  0.63 & (0.56, 0.73)\\
\hline
\end{tabular}
\end{center}
\end{table}





\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Bootstrap_PDF_H12T12.pdf}
\end{center}
\caption{Empirical density of the predicted 90\% fractile value
based on sample of size 20 for the sea state with $H_s=12$~m,
$T_p=12$~s. The $*$ indicates the limits of CI$_{0.95}$.}
\label{fig:11.7.f2e}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Bootstrap_PDF_H14_7T16_5.pdf}
\end{center}
\caption{Empirical density of the predicted 90\% fractile value
based on sample of size 20 for the sea state with $H_s=14.7$~m,
$T_p=16.5$~s. The $*$ indicates the limits of CI$_{0.95}$.}
\label{fig:11.7.f2f}
\end{figure}




Let us now compare the results provided by the Gumbel method\index{Gumbel method}
previously discussed, with the results obtained by the point process
method.\index{point process method} Using a Levenberg-Marquardt least squares
optimization method, cf. Chapter~\ref{chap:acer}, leads to the results shown in Figs.~\ref{fig:11.7.f42} and \ref{fig:11.7.f52}.




\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Kvitb.pdf}%\center
\end{center}
\caption{Log plot of the empirical and fitted upcrossing rate with
the reanchored 95\% empirical confidence band ($- -$) and fitted
confidence band (-- $\cdot$) for the sea state with $H_s=12$~m,
$T_p=12$~s. 90\% fractile estimate = 0.47~m, with CI$_{0.95}$ =
(0.45, 0.50).} \label{fig:11.7.f42} 
\end{figure}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=3in]{Figures/Kvitb3.pdf}%\center
\end{center}
\caption{Log plot of the empirical and fitted upcrossing rate with
the reanchored 95\% empirical confidence band ($- -$) and fitted
confidence band (-- $\cdot$) for the sea state with $H_s=14.7$~m,
$T_p=16.5$~s. 90\% fractile estimate = 0.62~m, with 95\%
CI$_{0.95}$ = (0.58, 0.65).} \label{fig:11.7.f52}
\end{figure}


As shown in Figures~\ref{fig:f3} and \ref{fig:f5}, when the mean upcrossing
rate is plotted on a logarithmic scale, the tails are closely
linear. This means that the associated extreme value distribution
can be expected to be similar to a Gumbel distribution, which
would correspond to exactly linear tails, cf. Eq.~(\ref{eq:6.30}).
Although it may not be obvious that the data plotted in
Figures~\ref{fig:11.7.f2b} and \ref{fig:11.7.f2c} come from a
distribution very close to a Gumbel distribution, the approximately linear exponential decay of the crossing rate strongly supports this assumption. This is yet
another indication of the usefulness of the mean upcrossing rate
function.



\begin{figure}[!t]
\begin{center}
\includegraphics[width=2.75in]{Figures/crH12T12.pdf}
\end{center}
\caption{Mean upcrossing rate statistics along with 95\%
confidence bands ($- -$) for the sea state with $H_s=12$~m,
$T_p=12$~s, and response standard deviation $\sigma=0.047$~m. $*$ : Monte Carlo; $- \, - \, -$ :
linear fit.} \label{fig:f3}
\end{figure}



\begin{figure}[!bht]
\begin{center}
\includegraphics[width=2.75in]{Figures/crH14_7T16_5.pdf}
\end{center}
\caption{Mean upcrossing rate statistics along with 95\%
confidence bands ($- -$) for the sea state with $H_s=14.7$~m,
$T_p=16.5$~s, and response standard deviation $\sigma=0.068$~m. $*$ : Monte Carlo; $- \, - \, -$ :
linear fit.} \label{fig:f5}
\end{figure}





Aiming at $T=3$-hour extreme response prediction, one needs
upcrossing rates down to about $10^{-4} - 10^{-6}$. Accurate
estimates based on direct Monte Carlo simulation down to this
order are expensive in terms of CPU time for the considered structure. It is therefore
convenient when accurately estimated upcrossing rates down to
about $10^{-3}$ can be used as a basis for extrapolation down to
the appropriate response level $\xi$ (with
$\nu^+(\xi)\approx10^{-5}$), as illustrated in
Figures~\ref{fig:11.7.f42} - \ref{fig:11.7.f52}.

Returning now to the specific prediction of the 90\% fractile of
the 3-hour extreme value distribution, $L_{CR}$,
Figures~\ref{fig:11.7.f42} and \ref{fig:11.7.f52} lead to the
estimates listed in Table~\ref{tab:11.7.3}. The estimated 95\%
confidence intervals are given in Table~\ref{tab:11.7.3}, and indicated in the figures. They are significantly
smaller than those obtained by the Gumbel method. The prediction accuracy is thus
significantly higher for the proposed method. However, it is also
observed that there is good agreement between the $L_{CR}$-values and the $L_G$-values, which is to be expected because the exact extreme value distribution is very close to a Gumbel distribution.

\bigskip

\begin{table}[!hbt]
	\begin{center}
		\caption{90\% fractile values $L_{CR}$ of the fitted Gumbel distributions }\label{tab:11.7.3}
		\begin{tabular}{c c c c}
			\hline
			$H_s$ (m) & $T_p$ (s) &  $L_{CR}$ (m) & CI$_{0.95}$ \\
			\hline
			12.0 & 12.0 & 0.47 & (0.45, 0.50)\\
			14.7 & 16.5 &  0.62 & (0.58, 0.65)\\
			\hline
		\end{tabular}
	\end{center}
\end{table}




\subsection{Combination of multiple stochastic load effects}\label{ssec:slc}\index{load effects!combination}\index{combination of load effects}


A prominent problem in the design of structures subjected to
random loads is to find methods for the combination of resulting
load effects at high and extreme response levels. In codified
design, this is usually implemented as linear combination rules of
specified characteristic values of the individual load effects
\citep{makrli86,melc99}. For nonlinear dynamic structures, the
precision level of such procedures would seem highly questionable.
One of the reasons for adopting such simplified procedures is the
complexity of the task to accurately predict the extreme value
statistics of the combined load effects, even in the case of
linear combinations. Over the years, several simplified procedures
have been suggested for the linear combination of load effects,
most notably the Ferry Borges-Castanheta method\index{Ferry Borges-Castanheta method} \citep{dima96},
Turkstra's rule\index{Turkstra's rule} \citep{turk70,makrli86}, the load coincidence
method\index{load coincidence method} \citep{wen90,melc99}, the SRSS method\index{square root of sum of squares
(SRSS) method} \citep{wepe81,wen90}, and the point crossing
approximation method\index{point crossing
approximation method} \citep{laco81,makrli86}. A main shortcoming
of these combination procedures is that they apply only to the
case of independent load effect components. A method for lifting
this restriction from the point process approximation is proposed
by \cite{toro84}. An effort to extend Turkstra's rule to
dependent processes is described by \cite{naro00}. In this
section, the use of the point process method for
stochastic load effect combination problems is illustrated.


The general formulation of the load effect combination problem to
be studied here is the following:
\begin{equation}\label{eq:H}
H(t)=h[X_1(t), \ldots, X_N(t)],
\end{equation}
where the stochastic load effect component processes $X_1(t),
\ldots, X_N(t)$ are combined according to a specified
deterministic function $h$ to produce the load effect combination
process $H(t)$. The component processes may derive from a
vector solution process of a dynamic model for the structural
response of an offshore platform to random waves. They may often
be modeled as stationary stochastic processes, but that is not a
requirement for the application of the point process method.

The typical problem to answer concerning the load effect
combination process $H(t)$ is to determine the probability of
exceeding a critical threshold $h_c$ during a specified time
interval $T$. Let us call this the failure probability and denote
it by $p_f = p_f(T)$. Defining $M(T) = \max\{H(t): 0 \le t \le T
\}$, it is realized that the goal is to calculate,
\begin{equation}\label{eq:H0}
p_f = 1 - \Prob( M(T) \le h_f ) \, .
\end{equation}

In many practical applications, the structure of the process
$H(t)$ is quite involved and the dimension $N$ can be high. This
makes a direct analytical approach virtually impossible. In such cases, Monte Carlo simulations of some sort would
seem to be the most attractive way to provide estimates of the
failure probability.


Two different load effect combination examples are used for
illustration purposes: von Mises yielding stress\index{von Mises yielding stress} and linear
combination of non-Gaussian load effects.\index{combination of load effects!non-Gaussian} In both cases, the load
components are correlated stochastic processes. A Newmark
integration method was used to produce accurate response time
series, cf. \cite{namo13}.

The load effect components $X_i(t)$ are modeled as stationary
processes, being the response of Duffing-type systems to the same
stationary Gaussian white noise excitation $W(t)$ with
$\EE[W(t)W(t+\tau)] = \delta(\tau)$, where $\delta(\cdot)$ denotes the Dirac delta function. That is,
\begin{equation}\label{eq:Dynamic}
\ddot{X}_i + 2\zeta_i\omega_i\dot{X}_i + \omega_i^2(X_i +
\varepsilon X_i^3) = W(t)/m_i \,,
\end{equation}
with specific damping constants $\zeta_i$ and resonance
frequencies $\omega_i = 2\pi/T_i$, and $m_i$ represent masses,
$i=1, \ldots, N$. In this subsection $N$ is taken to be 3.

First, the linear case, i.e., $\varepsilon=0$ in
Eq.~(\ref{eq:Dynamic}), is considered. In this case, the (two-sided) PSD of the
process $X_i(t)$ is given as $S_i(\omega) = |A_i(\omega)|^2$,
where,
\begin{equation}\label{eq:PSD}
A_i(\omega) = \frac{1}{ \sqrt{2 \pi} m_i
(-\omega^2+2 \ii \zeta_i\omega_i\omega+\omega_i^2) }\,,
\end{equation}
with $i=1,2,3$. The correlation coefficient
$\rho_{ij}=E[X_i(t)X_j(t)]/\sigma_i\sigma_j$ between $X_i(t)$ and
$X_j(t)$ is given by,
\begin{equation}\label{eq:rho}
\rho_{ij} = \int_{- \infty}^\infty A_i(\omega)\, A_j^*(\omega) \,
d\omega / \sigma_i\sigma_j\,.
\end{equation}
Here, $\sigma_i^2=\textrm{Var}[X_i(t)]=\int_{- \infty}^\infty
S_i(\omega)d\omega$ is the variance of $X_i(t)$, $i=1, 2, 3$.



\subsubsection*{von Mises stress combination}\label{num1}\index{von Mises stress combination}

Let $\sigma_x, \sigma_y, \sigma_z$ be axial stresses and
$\tau_{xy}, \tau_{xz}, \tau_{yz}$ be shear stresses in a
structural element. According to the von Mises yield criterion
\citep{makrli86}, yielding occurs if,
\begin{equation}\label{eq:Mises0}
(\sigma_x-\sigma_y)^2+(\sigma_x-\sigma_z)^2+(\sigma_y-\sigma_z)^2
+6(\tau_{xy}^2+\tau_{xz}^2+\tau_{yz}^2)\geq2\sigma_Y^2 \,,
\end{equation}
where $\sigma_Y$ is the yield stress. In many cases, in practice,
several stress components are zero. In the analysis carried out
here, only $xy$-plane stresses are encountered, meaning that
$\sigma_z=\tau_{xz}=\tau_{yz}=0$. According to the von Mises
criterion, yielding then occurs if,
\begin{equation}\label{eq:Mises}
(\sigma_x-\sigma_y)^2+\sigma_x^2+\sigma_y^2+6 \tau_{xy}^2 \geq
2\sigma_Y^2.
\end{equation}

The load effect vector in three-dimensional space is introduced as,
\begin{equation}\label{eq:Mises1}
(X_1(t), X_2(t), X_3(t)) \propto (\sigma_x(t), \sigma_y(t) \,,
\tau_{xy}(t))\,,
\end{equation}
where the components $X_k(t)$ are determined by
Eq.~(\ref{eq:Dynamic}).


The von Mises criterion (Eq.~(\ref{eq:Mises})) states that yielding
(failure) occurs when,
\begin{equation}\label{eq:HMises}
H_{vM}(t) = \big( (X_1-X_2)^2+X_1^2+X_2^2+6X_3^2 \big)^{1/2} \ge
h_f \,,
\end{equation}
where $h_f \propto \sqrt{2}\sigma_Y$.



In the linear case, ($\varepsilon=0$ in
Eq.~(\ref{eq:Dynamic})), the process $H_{vM}(t)^2$ becomes a
quadratic expression in correlated stationary Gaussian processes. Hence, the saddle
point method described by \cite{naka04pem} is applicable. Because
this method is very accurate, it provides an opportunity to check
the accuracy and efficiency of the point process method.
Table~\ref{tab:T2} lists the damping values and resonance periods
used in Eq.~(\ref{eq:Dynamic}) for the load effect processes $X_i(t)$.

\begin{table}[!htb]
\begin{center}
\caption{Model parameters}\label{tab:T2}
\begin{tabular}{c c c c}
\hline	
$i$ & $\zeta_i$ & $T_i$ (s) & $m$ \\
\hline
1 & 0.04 & 1.8 & 1\\
2 & 0.04 & 2.0 & 1\\
3 & 0.04 & 2.2 & 1\\
\hline
\end{tabular}
\end{center}
\end{table}


For the linear case ($\varepsilon=0$), it is straightforward to
calculate the standard deviations and correlation coefficients for
the $X_k$, cf. Eq.~\eqref{eq:rho}: $\sigma_1=0.38$, $\sigma_2=0.45$, $\sigma_3=0.52$,
$\rho_{12}=0.36$, $\rho_{13}=0.14$, $\rho_{23}=0.41$. There is significant correlation between some of the load effect components. For each example, 20 time series of length
approximately 0.92 hours each were simulated. The total
computation time for each example was less than a minute,
including simulation time and optimization.


The log plot presented in Figure~\ref{fig:CR0} shows the optimized
fitted parametric curve for the case of linear dynamics, i.e.,
$\varepsilon = 0$, together with the confidence band generated by
the fitted parametric curves to the borders of the
reanchored empirical 95\% confidence band. For illustration
purposes, the predicted value given by
$\nu^+(\xi) = 10^{-6}$, which corresponds to the 99\% fractile
value of a 3-hour extreme value distribution is indicated. Predicted value =
$6.84 = 10.03\, \sigma$, with 95\% CI = (6.60, 7.12).



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_Mises_e=0.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{vM}(t)$:
Monte Carlo ($\bullet$), reanchored empirical 95\% confidence band
(- -), fitted curve (--~--), fitted confidence band (-- $\cdot$).
Linear case: $\varepsilon=0$. $a = 0.795$, $b = 0.996$, $c =
1.607$, $\ln q = - 0.249$, $\sigma=0.682$.} \label{fig:CR0}
\end{figure}

In Figure~\ref{fig:CR}, it is demonstrated that when the fitted
parametric curve shown in Figure~\ref{fig:CR0} is replotted as an
optimal transformed plot, which yields a straight line, the empirical curve is also largely indistinguishable from a straight line, verifying the validity of our assumption about the
representation of the upcrossing rate function. The results from
saddle point calculations, which are practically exact, are also 
plotted in Figure~\ref{fig:CR}, and they are seen to agree very
well with the extrapolated straight line results.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_Mises_e=0.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{vM}(t)$: Monte Carlo ($\bullet$), empirical confidence band (-
-), saddle point results (o), linear extrapolation (--~--). Linear
case: $\varepsilon=0$. $b = 0.996$, $\ln q = - 0.249$,
$\sigma=0.682$.} \label{fig:CR}
\end{figure}


Figure~\ref{fig:CRnl0} shows the optimized fitted parametric curve
for the case of nonlinear dynamics, i.e., $\varepsilon = 1$,
together with the confidence band as generated previously. Predicted
99\% fractile value of 3-hour extreme = 4.29 = 8.04 $\sigma$,
95\% confidence interval = (4.21, 4.36). Figure~\ref{fig:CRnl}
presents the mean upcrossing rate function of $H_{vM}(t)$ plotted
on the optimal transformed scale. It is observed that the
assumption is again fully verified.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_Mises_e=1.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{vM}(t)$:
Monte Carlo ($\bullet$), reanchored empirical 95\% confidence band
(- -), fitted curve (--~--), fitted confidence band (-- $\cdot$).
Nonlinear case: $\varepsilon=1.0$. $a = 0.082$, $b = -0.382$, $c =
3.334$, $\ln q = 0.124$, $\sigma=0.534$.} \label{fig:CRnl0}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_Mises_e=1.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{vM}(t)$: Monte Carlo ($\bullet$), empirical confidence band (-
-), linear extrapolation (--). Nonlinear case: $\varepsilon=1.0$.
$b = -0.382$, $\ln q = 0.124$, $\sigma=0.534$.} \label{fig:CRnl}
\end{figure}




\subsubsection*{Linear combination of non-Gaussian load effects}\label{num2} \index{combination of load effects!non-Gaussian}

To get a flexible model that also provides a convenient way of
investigating the effect of statistical dependence between load
components, an example from \cite{naro00} is used. The input
load components are again assumed to be given by
Eq.~(\ref{eq:Dynamic}), with the same parameters as in the
previous case study, cf. Table~\ref{tab:T2}.


As an example of non-Gaussian load effect component processes, memoryless transformations of the input processes $X_i(t)$ provided by Eq.~(\ref{eq:Dynamic}) is considered. In this section, it is
assumed accordingly that,
\begin{equation}\label{eq:Z}
H_{lc}(t)= \sum_{i=1}^3 Z_i(t),
\end{equation}
with
\begin{equation}\label{eq:Zj}
Z_i(t)= X_i(t)|X_i(t)|^{\alpha}, \mbox{ } -1<\alpha<1.
\end{equation}
Two $\alpha$-values were chosen here, $- 0.5$ and $0.5$. The number
of terms in the sum (\ref{eq:Z}) is chosen to be three, but it can be
any positive integer because it does not matter much for the Monte
Carlo simulation.

Figure~\ref{fig:CR01} shows the optimized fitted parametric curve
for $\alpha = - 0.5$ and the case of linear dynamics, i.e.
$\varepsilon = 0$, together with the confidence band generated by
the allowed parametric curves. The predicted 99\% fractile value of 3
hour extreme = 3.42 = 3.34 $\sigma$, with the 95\% confidence interval =
(3.36, 3.49).  In Figure~\ref{fig:CR1} it is demonstrated that
when the fitted parametric curve shown in Figure~\ref{fig:CR01} is
replotted as an optimal transformed plot, the empirical curve is largely indistinguishable from a
straight line, supporting the assumption about the representation of the
upcrossing rate function.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_LoadComb_e=0_p=n_5.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{lc}(t)$:
Monte Carlo ($\bullet$), empirical confidence band (- -), fitted
curve (-- --), fitted confidence band (-- $\cdot$). Linear case:
$\varepsilon=0$. $a = 0.147$, $b = 0.129$, $c = 3.742$, $\ln q =
-1.103$, $\sigma=1.026$, $\alpha=-0.5$.} \label{fig:CR01}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_LoadComb_e=0_p=n_5.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{lc}(t)$: Monte Carlo ($\bullet$), confidence band (- -),
linear extrapolation (--). Linear case: $\varepsilon=0$. $b =
0.129$, $\ln q = -1.103$, $\sigma=1.026$, $\alpha=-0.5$.}
\label{fig:CR1}
\end{figure}


Figure~\ref{fig:CR02} shows the optimized fitted parametric curve
for $\alpha = 0.5$ and the case of linear dynamics, together with
the confidence band generated by the allowed parametric curves.
The predicted 99\% fractile value of 3 hour extreme = 5.15 = 7.69
$\sigma$, with the 95\% confidence interval = (4.95, 5.37). In
Figure~\ref{fig:CR2} it is demonstrated that when the fitted
parametric curve shown in Figure~\ref{fig:CR02} is replotted as an
optimal transformed plot, the empirical curve is largely indistinguishable from a
straight line, again supporting the assumption about the representation of the
upcrossing rate function.



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_LoadComb_e=0_p=p_5.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{lc}(t)$:
Monte Carlo ($\bullet$), empirical confidence band (- -), fitted
curve (-- --), fitted confidence band (-- $\cdot$). Linear case:
$\varepsilon=0$. $a = 1.771$, $b = 0.449$, $c = 1.280$, $\ln q =
-0.980$, $\sigma=0.670$, $\alpha=0.5$.} \label{fig:CR02}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_LoadComb_e=0_p=p_5.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{lc}(t)$: Monte Carlo ($\bullet$), confidence band (- -),
linear extrapolation (--). Linear case: $\varepsilon=0$. $b =
0.449$, $c = 1.280$, $\ln q = -0.980$, $\sigma=0.670$,
$\alpha=0.5$.} \label{fig:CR2}
\end{figure}



Figures~\ref{fig:CR03} and \ref{fig:CR3} show the log plot
and the optimal transformed plots for the case $\alpha = -0.5$,
under nonlinear dynamics with $\varepsilon = 1$. The predicted 99\%
fractile value of the 3-hour extreme = 3.04 = 3.18 $\sigma$, with the 95\%
confidence interval = (3.00, 3.09). By way of comment to the
rather high optimal value for $c$ that was found in this
particular case, it may be of interest to observe that the
function $F(q^*(c), a^*(c), b^*(c), c)$ is almost constant for a
range of values from about $c = 4$ to 7, where the
calculations stopped. From the transformed optimal plot, the assumption is also still verified for this case.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_LoadComb_e=1_p=n_5.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{lc}(t)$:
Monte Carlo ($\bullet$), empirical confidence band (- -), fitted
curve (-- --), fitted confidence band (-- $\cdot$). Nonlinear
case: $\varepsilon = 1$. $a = 0.014$, $b = -0.456$, $c = 5.423$,
$\ln q = - 1.048$, $\sigma = 0.958$, $\alpha=-0.5$.}
\label{fig:CR03}
\end{figure}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_LoadComb_e=1_p=n_5.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{lc}(t)$: Monte Carlo ($\bullet$), confidence band (- -),
linear extrapolation (--). Nonlinear case: $\varepsilon = 1$. $b
= -0.456$, $\ln q = - 1.048$, $\sigma = 0.958$, $\alpha=-0.5$.}
\label{fig:CR3}
\end{figure}


Figures~\ref{fig:CR04} and \ref{fig:CR4} show the log plot
and the optimal transformed plots for the case $\alpha = 0.5$,
under nonlinear dynamics with $\varepsilon = 1$. The predicted 99\%
fractile value of 3-hour extreme = 3.06 = 6.147 $\sigma$, with the 95\%
confidence interval = (2.97, 3.12). From the transformed optimal
plot, it is again seen that the assumption is fully verified also
for this case.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/conf_LoadComb_e=1_p=p_5.pdf}%\center
\end{center}
\caption{Log plot of the mean upcrossing rate of $H_{lc}(t)$:
Monte Carlo ($\bullet$); empirical confidence band (- -); fitted
curve (-- --); fitted confidence band (-- $\cdot$). Nonlinear
case: $\varepsilon = 1$. $a = 1.872$, $b = 0.030$, $c = 1.767$,
$\ln q = - 0.527$. $\sigma = 0.498$, $\alpha=0.5$.}
\label{fig:CR04}
\end{figure}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/fit_LoadComb_e=1_p=p_5.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the mean upcrossing rate of
$H_{lc}(t)$: Monte Carlo ($\bullet$), confidence band (- -),
linear extrapolation (--). Nonlinear case: $\varepsilon = 1$. $b
= 0.030$, $\ln q = - 0.527$. $\sigma = 0.498$, $\alpha=0.5$.}
\label{fig:CR4}
\end{figure}




20 maxima $M_j, j=1,..,20$ are extracted, one from each
realization, in order to view it on a Gumbel plot. The latter is a
plot of $M_j$ versus $-\ln(\ln\frac{20+1}{j})$, $j=1,..,20$.

A 95\% confidence interval for the response
level $L_{90}$ of the Gumbel distribution not being exceeded during time
$T=500\max(T_1,T_2,T_3)$ (cf. Table~\ref{tab:T2}), with probability 90\% based on a sample of
size 20 can be obtained by the MC technique and with
parametric bootstrapping from the fitted Gumbel distribution
\citep{dahi97}. One hundred thousand bootstrap samples where used to estimate
the density of the 90\% fractile, and from this density, the desired
confidence interval was extracted. Figures~\ref{fig:CDFnl} and
\ref{fig:Gnl} present the $L_{90}$ estimates by the point process
and Gumbel methods for the nonlinear system $\varepsilon=1$,
$\alpha=0.5$. Figure~\ref{fig:Bnl} presents the parametrically
bootstrapped density of $L_{90}$ from the Gumbel distribution.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/MCnl0_5.pdf}%\center
\end{center}
\caption{Distribution by the point process method; nonlinear system,
$\varepsilon=1$, $\alpha=0.5$.} \label{fig:CDFnl} 
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/Gnl0_5.pdf}%\center
\end{center}
\caption{Gumbel plot of 20 maxima; nonlinear system,
$\varepsilon=1$, $\alpha=0.5$.} \label{fig:Gnl} 
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/BootstrNl0_5.pdf}%\center
\end{center}
\caption{Bootstrapped density; nonlinear system, $\varepsilon=1$,
$\alpha=0.5$.} \label{fig:Bnl} 
\end{figure}


The estimate of $L_{90}$ from the point process method is 3.70 with
95\% confidence interval $(3.65, 3.75)$. For the Gumbel method, the
estimate of $L_{90}$ is 3.55 with 95\% confidence interval $(3.23,
3.85)$.





\subsection{Total surge response of a TLP}\label{ssec:tsrtlp} \index{surge response of a TLP}

The next example illustrates the problem of predicting the total surge
response of a tension leg platform (TLP) in random waves. The TLP concept was developed for production of oil at offshore fields. A simple rendition of a TLP structure is presented in Fig.~\ref{fig:tlp}. 
With the tools developed in this
chapter, the problem of response prediction may be solved in a rather
satisfactory manner. The presented material is largely taken from
\cite{nagate07aor}. A Monte Carlo based approach to the investigation of the response statistics of an offshore structure is also presented in \cite{saganali11oe}.

\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.75in]{Figures/tlp2.pdf}%\center
	\end{center}
	\caption{A sketch of a TLP structure} \label{fig:tlp} 
\end{figure}


The equations of motion for a floating, rigid-body TLP structure
subjected to environmental forces such as wind, waves and current
would generally be written as,
\begin{equation}\label{eq:dm}
\pmb{M} \, \ddot{\pmb{Z}}(t) + \pmb{H}( \pmb{Z}(t),
\dot{\pmb{Z}}(t), t) = \pmb{F}(t)\,.
\end{equation}
Here, $\pmb{M}$ denotes a generalized $6 \times 6$ mass matrix,
$\pmb{Z} = (Z_1, \ldots, Z_6)^T$ = the structure's response
vector, while $\pmb{H}$ denotes a nonlinear vector function.
$\pmb{F}(t)$ denotes a stochastic loading process, which in
general also depends on the response of the structure. This point
is further discussed later.

Because the main purpose of this section is to illustrate the
versatility and accuracy of the proposed method, we chose to
discuss a simplified SDOF model for the
surge response of the TLP in random waves. Except for interaction
effects between different motion modes, the SDOF model allows for
the introduction of most of the relevant nonlinear effects that
may influence the response of the TLP. Hence, the following SDOF
equation of motion is studied,
\begin{equation}\label{eq:dm2}
{M} \, \ddot{{Z}}(t) + {H}( {Z}(t), \dot{{Z}}(t), t) = {F}(t),
\end{equation}
where $Z=Z(t)$ denotes the surge response of the TLP; $M$ is the
mass of the platform, including added mass; and $H$ is a nonlinear
function to be specified. As discussed in Chapter~\ref{chap:ch},
the hydrodynamic loading process $F(t)$ is assumed to consist
of two components: a linear, first-order (wave frequency) term
$F_1(t)$, and a nonlinear, slowly varying, second-order term
$F_2(t)$.

To set up the proper dynamic model for the surge response
$Z(t)$, it is necessary to take into account the fact that
hydrodynamic loading on a floating body depends on the motions of
the body. In case of the slow drift motions of the TLP, it is of some importance to take into account the
dependence of the slow drift force $F_2(t)$ on the slowly varying
surge velocity $\dot{Z}_2(t)$. Because of a nonlinear dynamic
model, a definition of the slow drift
response $Z_2(t)$ has to introduced. A suitable definition would seem to be the
following: The slow-drift response $Z_2(t)$ is obtained from the
total response $Z(t)$ by a low-pass filter that removes all wave
frequency components. In practice, this may be achieved by a
running mean operator used iteratively.

To account for the dependence of the slow drift force on the
slowly varying velocity, it is appropriate to write
$F_2(t,\dot{Z}_2(t))$ rather than $F_2(t)$. However,
$F_2(t,\dot{Z}_2(t))$ is not directly available to us, but only
$F_2(t)\equiv F_2(t,0)$ given by an equation entirely similar to Eq.~(\ref{eq:cf7}). Since, in the
context of slow drift motions, $\dot{Z}_2(t)$ is small, the
following approximation is adopted,
\begin{equation}\label{eq:slow drift damping}
F_2(t,\dot{Z}_2(t)) \approx  F_2(t,0) + \frac{\partial
F_2(t,0)}{\partial \dot{Z}_2} \, \dot{Z}_2(t).
\end{equation}
It is shown by \cite{najo93pem} that for a TLP structure
$\frac{\partial F_2(t,0)}{\partial \dot{Z}_2} \approx -c \, F_2(t,
0)\equiv - c \, F_2(t)$ for a suitable constant $c>0$ may to some
extent serve as a useful approximation to capture qualitatively
the time-variant damping effect, which is the result of the
expansion in Eq.~(\ref{eq:slow drift damping}).

The first dynamic model adopted for $Z(t)$ is now the following,
\begin{equation}\label{eq:ds}
\tilde{M} \, \ddot{Z}+C \,\dot{Z} + K (Z+\varepsilon Z^3) = F_1(t)
+ F_2(t,\dot{Z}_2) \approx F(t) - cF_2(t)\, \dot{Z}_2.
\end{equation}
Here, $\tilde{M} = M + \tilde{m}$, where $\tilde{m}$ is an
appropriately chosen (constant) added mass. $C$, $K$, and
$\varepsilon$ are suitably chosen positive constants. This
equation is rewritten in the equivalent form,
\begin{equation}\label{eq:ds1}
\ddot{Z}+2\omega_e \, \zeta \dot{Z} + 2\omega_e \tilde{c} F_2(t)
\, \dot{Z}_2 + \omega_e^2(Z + {\varepsilon}Z^3) =
\frac{F(t)}{\tilde{M}},
\end{equation}
where $\omega_e^2 = K/\tilde{M}$, $\zeta = C/(2 \omega_e
\tilde{M})$, and $\tilde{c} = c/(2 \omega_e \tilde{M})$.

Thus, the dynamic system is nonlinear with a Duffing-type
hardening stiffness nonlinearity and time-varying damping. 
For the TLP the relative damping coefficient $\zeta$ is usually
small. As a consequence, the contribution from the time-varying term is
non-negligible, especially for severe seas for which the slow
drift response is significant. The third-order term in the
restoring force, generally referred to as the \emph{set down} effect,
is caused by the fact that the tethers will induce the TLP to act
like an inverted pendulum. Note that the set down effect will also have
an influence on the hydrodynamic loading process, which
depends not only on $\dot{Z}$, but also on $Z$. Even
if this dependence could have been taken into account, it was neglected here because it is of minor importance.


For the numerical simulations, a particular model of a TLP is
considered, and the corresponding LTF and QTF are computed using the second-order diffraction program \citep{wamit}. For simplicity, unidirectional
seas are used, meaning that the directional argument $\beta$ is
skipped. This simplification should have no effect on the
conclusions based on the comparison of accuracy. The combined
first-order and second-order slowly varying surge deck motion is
studied applying the single-degree-of-freedom model. The TLP
particulars are detailed in Table~\ref{tab:1}, and the subsurface part of the
structure is shown in Figure~\ref{fig:TLP}, but without the vertical tethers.

\begin{figure}[!hbt] %[yl,y]=TLPcopy(0)
\begin{center}
\includegraphics[width=2.75in]{Figures/mini_TLP_906.pdf}%\center
\end{center}
\caption{Sketch of the submerged part of the TLP. Units in meters}
\label{fig:TLP} 
\end{figure}

The values in Table~\ref{tab:1} are used to obtain the second-order response. This means that for the second-order response, a simplified version of Eq.~(\ref{eq:ds}) was used, where
mass $\tilde{M}$, stiffness $K$, and damping coefficient $C$ are
frequency independent, which is a good approximation for the slow
drift motion. The time-invariant damping part $\zeta$ is
considered to be $5$\%.


\begin{table}[!htb]
\centering
\caption{Particulars of the TLP}\label{tab:1}
\begin{tabular}{l r}
\hline
Column diameter $D$ (m) & 10.0\\
Eigenperiod surge $T_e$ (s) &  128.8\\
Relative damping $\zeta$  & 0.05 \\
Total mass (incl. added mass) $\tilde{M}$ (kg) & $1.5\cdot10^7$\\
\hline
\end{tabular}
\end{table}


Two versions of Eq.~(\ref{eq:ds1}) are used. The first version
is a linear, time-invariant model obtained by putting $\tilde{c} =
\varepsilon = 0$. The second version is the fully nonlinear model
where the parameter $\tilde{c}$ in Eq.~(\ref{eq:ds1}) is chosen such
that $\Var[2\omega_e \tilde{c} F_2(t) \, \dot{Z}_2(t)]$ is about
10\% of $\Var[2\omega_e \, \zeta \dot{Z}(t)]$.  The parameter
${\varepsilon}$ is estimated from the condition that $0.2 Z(t)\geq
{\varepsilon} Z^3(t)$ when $Z(t)\leq 6\sigma_Z$, i.e. even in the
extreme response region, stiffness hardening contributes not more
than 20\% relative to the linear part for severe sea, which lead
to ${\varepsilon}=1.36 \cdot 10^{-4}$. Finally, the following
approximate values were found: $\tilde{c}=30/(\tilde{M}g)$ for moderate seas, and
$\tilde{c}=90/(\tilde{M}g)$ for severe seas, where $g=9.81\textrm{m}/\textrm{s}^2$, cf. Eq.~\eqref{ex2}. The adopted parameter values are largely arbitrary, but the choices made seem to provide
a reasonable model for the chosen TLP structure.

To get an accurate representation of the response process, there
is a specific requirement that must be observed. Because the damping
ratio is only 5\%, the frequency resolution $\Delta\omega$ must
secure a sufficient number of frequency values over the resonance
peak. This will ensure that the second-order, difference-frequency
response component captures the TLP surge dynamics with sufficient
accuracy. It is commonly required that there are at least 5
discrete frequencies over the frequency range where
$|\hat{L}(\omega)|^2$ is equal to or higher than half of the
resonance peak height $\max(|\hat{L}(\omega)|^2)$, where,
\begin{equation}
\hat{L}(\omega) = \big(-\omega^2 + 2 \ii \zeta \,\omega_e \,\omega
+ \omega_e^2\big)^{-1}.
\end{equation}
For the surge force QTF $\hat{K}_2(\omega,\omega')$, a suitable
initial frequency grid must be chosen for which the values of the
force QTF are calculated. The calculation of the force QTF is
generally the most time-consuming part of the numerical
analysis. Therefore, the initial grid is usually rather coarse to
avoid excessive computer time. For the calculations at hand, the
discrete frequency range was the following: $\omega_1 = 2\pi/30.0,
\ldots, \omega_n = 2\pi/4.0$ (rad/s), $n=30$. This necessitates
the use of an interpolation procedure to be able to provide values
of the QTF on a much finer grid than the initial one to comply
with the requirement of sufficient frequency resolution to capture
the dynamics of slow drift motion. In this chapter cubic spline
interpolation is used. In the particular case considered
here, the resolution requirement led to the choice
$\Delta\omega=0.0018$ rad/s and $L=760$ interpolated discrete
frequencies.

The random stationary sea state is specified by a JONSWAP
spectrum, which is given as follows,
\begin{equation}\label{ex2}
S_\eta(\omega)=\frac{\alpha g^2}{\omega^5} \exp\Bigg\{ -
\frac{5}{4} \Big(\frac{\omega_p}{\omega}\Big)^4 + \ln \kappa \,
\exp\Big[ - \frac{1}{2 \chi^2} \Big(\frac{\omega}{\omega_p} - 1
\Big)^2 \Big] \Bigg\},
\end{equation}
where $g = 9.81$ ms$^{-2}$, $\omega_p$ denotes the peak frequency
in rad/s, and $\kappa$ and $\chi$ are parameters affecting the spectral shape.
 $\chi = 0.07$ when $\omega \le \omega_p$, and $\chi = 0.09$ when
 $\omega > \omega_p$. The parameter $\kappa$ is chosen to be equal
 to 3.3, which is a rather typical value. The parameter $\alpha$ is determined from the following
 empirical relationship,
\begin{equation}\label{ex3}
\alpha = 5.06 \Big(\frac{H_s}{T_p^2}\Big)^{2} \big( 1 - 0.287 \,
\ln \kappa \, \big)\,,
\end{equation}
where $H_s$ denotes the significant wave height and $T_p= 2 \pi/
\omega_p$ is the spectral peak wave period. Table~\ref{tab:16.6} presents the sea state parameters,
along with the corresponding response standard deviations. For more results, cf. \cite{nagate07aor}.


%C:\Comparison\TLPnonlin, MorisonNew
%for plotting plotTLP, plotTLP2

\begin{table}[!htb]
\centering
\caption{Representative sea state, along with response standard deviations for the linear and nonlinear TLP model.}\label{tab:16.6}
\begin{tabular}{c c c c}
\hline
$H_s$ (m) & $T_p$ (s) & $\sigma_Z$ (m), lin. & $\sigma_Z$ (m), nonlin. \\
\hline 10.0 & 11 & 9.3 & 8.2 \\
\hline
\end{tabular}
\end{table}
%\vspace{5mm}



Figure~\ref{fig:absLTF} shows the LTF for the wave exciting force
amplitude, while Figure~\ref{fig:absQTF} depicts the spline
interpolated QTF. Because the QTF is complex valued, only its
absolute value is plotted.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/LTF_r.pdf}%\center
\end{center}
\caption{Wave exciting force amplitude, surge LTF.}
\label{fig:absLTF} 
\end{figure}


\begin{figure}[!hbt] %[yl,y]=TLPcopy(0)
\begin{center}
\dofigs{3.4in}{2in}{Figures/QTF2.pdf}%\center
\end{center}
\caption{Wave exciting force amplitude, surge QTF.}
\label{fig:absQTF}
\end{figure}


For the chosen sea state, Figures~\ref{fig:CR1_1000} and
\ref{fig:CR1_25} present the corresponding response tail crossing
rates obtained by Monte Carlo simulation for the linear system
given by putting $\varepsilon = 0$ and $\tilde{c}  = 0$ in
Eq.~(\ref{eq:ds1}). Figure~\ref{fig:CR1_1000} shows the results obtained from 1,000
realizations, which requiring less than 1~hour on a laptop computer. The crossing rate plots are
done on the transformed scale, see Eq.~(\ref{eq:11.3.8}). Extreme
response prediction based on Eq.~(\ref{eq:11.3.1}) will typically
involve crossing rates of the orders $10^{-6}-10^{-7}$, but to
illustrate the achieved accuracy the extrapolated results at the
crossing rate level $10^{-10}$ are highlighted.
Figure~\ref{fig:CR1_1000} also shows the
highly accurate results obtained by using a saddle point
integration technique\index{saddle point integration technique} \citep{nakate06aor}. These results cannot be
distinguished from those obtained by linear extrapolation of the
mean upcrossing rate function provided by Monte Carlo simulations.
Hence, linear extrapolation can be done over several orders of
magnitude with high accuracy.

To illustrate the fact that good accuracy can be obtained with
much shorter time simulation records, Figure~\ref{fig:CR1_25} shows the results of using only 25~hours of simulated
response time histories, which required less than 10~minutes on a standard laptop PC.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/TLPlinM1000_bw.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the empirical crossing rates by Monte Carlo simulation ($*$)
with 95\% confidence bands ($--$) based on 1,000 hours of response time
histories for the case of linear dynamics ($\tilde{c}=
{\varepsilon}=0$). Saddle point integration results
($\frac{\qquad}{\qquad}$) coincide with the optimized linear fit
with $b=0.75\,\sigma_Z$, $q=0.0205$, where $\sigma_Z=9.3$~m (see
Table~\ref{tab:16.6}).  }
\label{fig:CR1_1000} \vspace{1cm}
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/TLPlinM25_bw.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the empirical crossing rates by Monte Carlo simulation ($*$)
with 95\% confidence bands ($--$) based on 25~hours of response time
histories for the case of linear dynamics ($\tilde{c}=
{\varepsilon}=0$); optimized linear fit ($\frac{\qquad}{\qquad}$)
with $b=0.75\,\sigma_Z$, $q=0.020$, where $\sigma_Z=9.3$~m (see
Table~\ref{tab:16.6}), for the case of linear dynamics ($\tilde{c}=
{\varepsilon}=0$).  }
\label{fig:CR1_25} \vspace{1cm}
\end{figure}




Figures~\ref{fig:CR10_1000} and \ref{fig:CR10_25} present response
tail crossing rates for the chosen sea state obtained by Monte Carlo
simulation for the fully nonlinear model given by
Eq.~(\ref{eq:ds1}). In this case the Monte Carlo simulation
results are the only results available for verification of the
extrapolation method for the nonlinear model, but the experience
from the linear case indicates that the results obtained from
Figure~\ref{fig:CR10_1000} are very
accurate. Again, the crossing rate plots are done on the
transformed scale, see Eq.~(\ref{eq:11.3.8}), and for illustration
purposes the extrapolated results at the crossing rate level
$10^{-10}$ are highlighted. These results can be compared with the
corresponding results obtained from only 25~h of simulated
response time histories shown in Figure~\ref{fig:CR10_25}. The agreement is again very good.



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/TLPnonlinM1000_bw.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the empirical crossing rates by Monte Carlo simulation ($*$)
with 95\% confidence bands ($--$) based on 1,000 hours of response time
histories for the fully nonlinear model; optimized linear fit
($\frac{\qquad}{\qquad}$) with $b=0.54\,\sigma_Z$, $q=0.0248$, where
$\sigma_Z = 8.2$~m (see Table~\ref{tab:16.6}). } \label{fig:CR10_1000} \vspace{1cm}
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.75in]{Figures/TLPnonlinM25_bw.pdf}%\center
\end{center}
\caption{Optimal transformed plot of the empirical crossing rates by Monte Carlo simulation ($*$)
with 95\% confidence bands ($--$) based on 25~hours of response time
histories for the fully nonlinear model; optimized linear fit
($\frac{\qquad}{\qquad}$) with $b=0.54\,\sigma_Z$, $q=0.0234$, where
$\sigma_Z = 8.2$~m (see Table~\ref{tab:16.6}). } \label{fig:CR10_25} \vspace{1cm}
\end{figure}




