Data-driven modeling of hydrodynamic loading using NARX and harmonic probing

Optimizing floating wind turbines and their mooring systems requires validated computational models that predict wave-frequency and low-frequency hydrodynamic loads. Low-frequency loads are crucial for determining extreme offsets and tension in mooring lines and are generally described by a quadratic transfer function. The quadratic transfer function, obtained with numerical tools, accurately predicts low-frequency loads in mild sea states. However, since the existing numerical methods are based on potential and perturbation theory, they generally fail to accurately predict low-frequency loads in moderate-to-extreme sea states where current, viscous, and beyond-second-order potential effects become significant. Developing a procedure for empirical transfer function estimation is, therefore, necessary to overcome these limitations. This paper describes an existing framework for estimating any higher-order transfer functions from experimental data. The framework employs a nonlinear auto-regressive model based on Kriging to establish a causal input/output relationship between the wave-elevation and hydrodynamic force histories exerted on the floater. Then, higher-order transfer functions are extracted using harmonic probing. The procedure was validated by estimating the linear surge transfer function of the INO WINDMOOR 12 MW floater using synthetic data. The data-driven results showed an excellent agreement with the theoretically computed transfer function.


Introduction 1.Background and Motivation
Accurate estimation of the wave loading on moored structures is essential for providing costcompetitive floating wind turbine designs.The challenge with these structures lies in their relatively soft mooring system which results in low natural frequencies in the horizontal degrees of freedom.Consequently, moored structures are sensitive to the low-frequency loads which contribute significantly to extreme offsets and tension in mooring lines.The low-frequency loads are generally described by a difference quadratic transfer function (QTF) which models the relationship between a bi-chromatic wave with frequencies f 1 and f 2 and the low frequency wave loading at |f 1 − f 2 |.Providing a reliable and computationally effective way of estimating the transfer function that models these loads is therefore of great importance for the offshore industry.However, as described in [1], computing an accurate estimation of the higherorder transfer functions is a challenging and computationally expensive task.Commercially available software such as Hydrostar [2] and WAMIT [3] do compute the first and secondorder transfer functions however they are based on potential flow and perturbation theory which inherently assume an inviscid, incompressible fluid and irrotational flow.It is important to remark that moderate-to-extreme sea states with steep waves violate the potential flow and perturbation theory assumptions which leads to inaccuracies in the transfer function (TF) estimation.Therefore, Morrison's equation with empirically calibrated coefficients is commonly implemented to account for the viscous effects that potential flow fails to capture.Alternatively, Computational Fluid Dynamics (CFD) models can be used to solve the wavestructure interaction problem.However their computational cost often makes them impractical in an industrial context where potential-flow-based solvers combined with Morrison's equation are still the preferred approach for solving this problem.A third alternative for TF estimation are data-driven methods.Although less known, data-driven approaches are not new and have been used for the purpose of transfer function estimation before.A prominent method is the cross-spectrum analysis with the cross-bi-spectral (CBS) method as an extension for estimating QTFs.Details about CBS are presented in [4], [5], and more recently in [1], while [6] applied CBS to the INO WINDMOOR floater and compared the results to potential theory.

Scope
An alternative approach within the data-driven domain are the nonlinear auto-regressive models with exogenous input (NARX).When used in tandem with harmonic probing, they allow for an estimation of any order transfer function.
The origin of NARX models dates back to the work of Billings in the 1980s when the ARMAX (auto-regressive moving average model with exogenous inputs) was extended to its nonlinear form i.e. the NARMAX model [7][8].This framework is the most generic and versatile formulation of a nonlinear discrete-time process that incorporates a noise model and is usually given a polynomial form.In [9] Billings showed several examples of how NARMAX naturally comes into existence when real continuous-time systems are discretized.A simplification to the NARMAX model is the NARX model which assumes that the noise process is white Gaussian [10].A neural-network-based NARX model was investigated in [11] together with harmonic probing.The study investigated different activation functions and validated the results against a nonlinear Duffing oscillator.Similarly, [10] used a Gaussian-process-based NARX i.e.Kriging-NARX, to extract an analytical expression of the higher-order frequency response functions of a Duffing oscillator.The obtained results showed a remarkable agreement.
The present paper shows the application of Kirging-NARX and harmonic probing framework for estimating the surge linear transfer function (LTF) of the INO WINDMOOR floater as a first step towards estimating the higher-order transfer function.The method employed here is based on the work presented by Worden et al. in [10][12] [13].To the authors' knowledge, this framework has not been validated for a moored floating structure before.It is worth noting that at present, the focus is cast solely on estimating the diffraction loading.Furthermore, it is assumed that the data related to the diffraction loading and wave elevation profile histories are available.Direct experimental measurements of diffraction loads in general is impractical.In contrast, the floater responses and wave elevation profile can be directly measured.Therefore, the diffraction loading is typically inferred via linear de-convolution of the floater response as presented by Sauder in [1].

Methodology
As mentioned earlier, the algorithm for estimating the transfer functions consists of two parts.Namely, the nonlinear auto-regressive model with exogenous input (NARX) and harmonic probing (HP), which will both be described in the following section.

Kriging-NARX
The main idea behind NARX (nonlinear auto-regressive model with exogenous input) is to construct and train a model (F) such that it can provide a prediction of the current output based on past input and output data. where Kriging-NARX is a special type of NARX that employs the concept of Kriging.Kriging dates back to the 1950s when it was first introduced for interpolating the profile of landscapes [14].Since then, the machine learning community has re-purposed Kriging for a variety of learning tasks.The book by Rasmussen and Williams [15] provides a good overview of this methodology, a consolidated summary is presented hereafter.
The underlying assumption of Kirging-NARX is that the model output is a realization of a Gaussian process (GP) [16].From the definition of a GP it follows that the distribution for a finite number of outputs follows a multivariate joint normal distribution.
(3) Notice that the mean of the joint distribution is assumed to be zero for simplicity.The method allows for a non-zero mean joint distribution though this adds additional model complexity and training parameters.In the above expression, the hyperparameter σ 2 e is the noise variance which in itself incorporates a nugget for numerical stability, whereas k(•, •) is a correlation function selected by the user.A number of correlation functions and their applications are listed in [17].A common choice is the squared exponential correlation function which is infinitely differentiable leading to smooth paths.This function was the natural choice for this study since the signals investigated here are of smooth nature.The squared exponential function was modified to include a scaling length on both input signals leading to the following expression: where ∥ ∥ is the Euclidean norm of the given vector difference while x f and x ζ refer to the autoregressive and exogenous part of the input vector, respectively.The expression above introduces three new hyperparameters: σ f and the scale lengths, θ f and θ ζ .Assume now that the training data is gathered as follows: with N representing the total number of data points.Then, the expression from (3) can be represented in a matrix form: The advantage of the joint Gaussian distribution is that the distribution of the unknown output f n for an unobserved input x n pre-conditioned on a set of observation {X, F} is also Gaussian.From ( 6), the corresponding mean and variance for the output f n have analytical expressions, which read, Notice above that the expectation is preconditioned on the past outputs F and also the training data X.That is, the training data become a part of the model equations and permanently "lives" in the algorithm.This feature provides an advantage of Kriging-NARX because it reduces the number of hyperparameters needed to capture the relationship of the modeled signals.In addition, this means that the model can be successfully trained with a relatively small amount of data.
The only remaining part is to close Eq.( 7), i.e. train the unknown hyperparameters, Training is performed by using a maximum likelihood estimation over a marginal evidence function.Put simply, since the outputs (or observations) F = {F(x 1 ), F(x 2 ), . . ., F(x N )} are assumed to have a joint normal distributed, maximizing their likelihood means finding the zero-mean joint normal distribution that best fits the data.For a given set of training data, (F, X), the likelihood function takes the following form: In practice, the hyperparameter values are obtained by minimizing the negative log likelihood function, which reads, Note that the dependence on the optimization parameters θ * is not explicitly shown above, but it is rather obvious from (4).Furthermore, the N 2 log(2π) term is left out since it does not affect the optimization.Due to the low number of hyperparameters to be optimized, the training can be completed with a simple gradient descent optimization algorithm.Note also from (10) that the optimization requires a matrix inversion.The matrix K(X, X) is of N × N size, where N is the number of observations in the training set.Since the numerical optimization algorithms rely on repeated evaluations of the objective function, the optimization may become costly for large data sets.Once the training is complete, the model (Eq.7) can be readily evaluated for a set of two new signals originating from the same underlying process for which the model was trained for.An application of this framework is presented in Section 3.

Harmonic Probing
Section 2.1 provided an efficient way to construct a nonlinear input/output relationship for a given system.However, in general, in vibration problems the transfer functions hold more information about the underlying physics of the problem i.e. which frequencies produce large outputs corresponding to resonance.Consequently, it is desirable to extract this physical quantity from the seemingly less informative NARX model.This is achieved with the method of harmonic probing (HP).HP was first introduced in [18] for continuous-time models.Since then, the method has been extended to discrete-time systems as well.A good overview is available in Chapter 8 of [12].A condensed description is presented hereafter.
The idea behind HP is to probe the system with simple harmonic excitations for which an analytical response is available.The analytical formulation between a set of simple harmonics and the corresponding response of the system is expressed through a Volterra series expansion.
where f 0 is a constant and: . . .
here f 1 is the response of the system, ζ is an input excitation, and h n is the n-th order impulse response function or also known as a Volterra kernel.The time-domain representation above can also be expressed in the frequency domain by simply taking the Fourier transform of Eqs.(11)(12)(13)(14).
where, F and Z are the frequency domain counterparts of the response and excitation while H n is the n-th order transfer function.Note that the expressions for F 2 . . .F n have been made symmetric at the expense of an additional integral, see [11].In theory, the Volterra expansion could capture the interactions of an infinite number of simple harmonics.In this study, the series was truncated after the second-order effects including only the linear H 1 (ω) and quadratic H 2 (ω 1 , ω 2 ) transfer functions.To obtain the final probing-friendly expression, the integrals in (17) can be evaluated analytically and the resulting expression transform back to the time domain.Similarly, the inverse Fourier transform is applied to ( 16) and together with (17) are substituted in (15).This gives the final time-domain response of the system which is obtained assuming a bichromatic excitation of amplitude unity: leading to the expression for the system response: where t n = n∆t.The time-step ∆t depends on the sampling of the NARX model.Determining the correct sampling rate depends on the dynamics of the modeled dynamic system.The equation above provides an analytical way to relate the response of the system to its transfer functions.Equation ( 20) can be equated to the linearized expectation of the predictor from (7).
This equality allows for the Kriging-NARX model from Section 2.1 to be related to a physicsbased representation i.e. the Volterra series expansion.The linearization above is necessary for converting the Kriging-NARX model into a polynomial form suitable for harmonic probing.
Note that the order of the transfer function of interest determines the minimum number of harmonics required as well as the number of terms included in the Taylor series expansion.Substituting the time-shifted input vector into the linearization of ( 21) gives a closed form equation where the only remaining unknowns are the transfer functions which can be isolated and solved for.The process of solving is completed in cascade starting from the lowest order and moving forward.First, the probed equation is solved for the constant terms leading to the calculation of the bias term H 0 .Then, a harmonic balance is solved for the frequency Ω 1 leading to the expression of the LTF, H 1 (Ω 1 ).Lastly, a harmonic balance is solved for the frequency Ω 1 + Ω 2 leading to the expression of the QTF, H 2 (Ω 1 , Ω 2 ).This recursive solving means that obtaining a p -th order transfer function requires the previous (p − 1) transfer functions to be solved for first since they appear in the coefficients of the harmonic balance.

Numerical Case Study
The framework described in the previous section was applied to establish a model of the surge wave loads, and limited for the time being to the surge linear transfer function (LTF), relating the wave elevation profile ζ to the hydrodynamic surge force f on the floater.A description of the floater under study is available in Section 3.1.The validation study used synthetic time-series data whose generation is described in Section 3.2.Lastly, the obtained results are presented in Section 3.4.

Floater Description
The floater used for this study was the INO WINDMOOR 12MW floater [19], a semi-submersible structure with three columns connected by pontoons.The structure has been modelled in Hydrostar [2] for the purposes of obtaining the transfer functions.The vessel mesh consisted of 11212 panels.Figure 1

Data Generation
The synthetic data was obtained by first producing a random wave elevation profile from a JONSWAP spectrum with a significant wave height of H s = 0.5 m and a period of T p = 13 s.A plot of the spectrum is shown in Figure 2. The surge transfer function of the floater, computed in  Hydrostar, together with a wave elevation profile, were used to generate a synthetic loading timeseries data using the approach outlined in Section 1.1.2 of [21].The Python toolbox Snoopy [22] was used for this purpose.The first-order force was calculated as follows: where |H(ω m )| and ϕ m are the amplitude and phase angle of the transfer function of the system, respectively.S(ω m ) is the spectral density of the JONSWAP spectrum, ∆ω is the frequency step and ϵ m is a random number between 0 and 2π.In practice, the loading time-series is more difficult to obtain and usually cannot be directly measured in the basin.One viable approach for obtaining the force signal from the measured response is outlined in [1] and [23].The main idea revolves around assimilating the measured displacement/velocity response history to a damped single-degree-of-freedom oscillator whose parameters have previously been estimated from a set of basin tests.The force history is then simply reconstructed by a linear de-convolution.

NARX Setting and Results
The synthetic time-series data for the wave elevation profile, ζ(t) and force, f (t) was divided into 20 segments where each segment consisted of 300 data points with a time-step, dt = s.Each segment was used for training a separate Kriging-NARX model.A crucial goal of the training process is to capture the underlying physics of the actual floater.In dynamics, modeling the behaviour of a structure, depending on its complexity, could require several past states.The INO WINDMOOR floater has a relatively complex wave-body interaction that requires a longer memory.The memory effects in the Kriging-NARX model are captured in the lagged input data.This study used n f = 20 and n ζ = 20 which corresponds to 24 s of lag on both the force and wave signals.It is worth noting that the wave loading process was assumed to be nearly causal.In reality, the wave excitation force and the wave elevation profile have a non-causal relationship however any effects of non-causality were not investigated in this study.The predictive capacity of each NARX was tested with a leave-one-out cross-validation technique.This process went on as follows, the 20 segments mentioned earlier were used as training data for the model producing 20 different NARX.Each one of the 20 NARX models was then validated against the remaining 19 segments.This ensured that each NARX was tested only on data that it has not seen before.The results for the 20 segments are displayed in the figure below.Figure 3 shows the output of 19 NARX models (blue) validated against a different (20th) reference data-set (blue).All models performed exceptionally well in estimating the force timeseries.An important thing to note here is that the predictive power of the model is limited to time-series sets corresponding to a particular sea state.For a new sea state, the model needs to be retrained again.

Transfer Function Results
Lastly, the harmonic probing algorithm was applied to the Kriging-NARX models resulting in a distinct transfer function for each model.The results are displayed in Figure 4.It is interesting to note that the region between 0.3−1.3rad/s, where most of the energy from the JONSWAP (see Figure 2) was located, exhibited lower uncertainty.A consequent analysis with a JONSWAP shifted towards the lower frequencies confirmed this observation.This explains the divergence of the results in the regions below 0.3 rad/s and above 1.5 rad/s.

Conclusion
This study applied and validated the Kriging-NARX and harmonic probing framework for extracting the diffraction load linear transfer function in surge of a floating structure.This method is tailored to be used with experimental data, although, this study utilized synthetic data generated for a given JONSWAP spectrum and a transfer function obtained with a potentialtheory-based software.For that reason, the obtained results reproduced the same transfer function that generated the loading time series.In reality, the obtained transfer function from experimental data is expected to differ than the one from potential-flow theory due to possible viscous effects that are present in the basin but are not captured with potential flow.The extent of this difference and overall validity of the model in severe sea states is yet to be further studied.For this sea state, the Kriging-NARX surrogate model provided excellent predictive performance for the specified training setting (i.e.number of lags and time step).Moreover, it is robust to overfitting when wave elevation profiles are all sampled from the same JONSWAP spectrum.The harmonic probing algorithm managed to successfully extract the linear transfer function from the trained wave-force relationship.It was observed that the uncertainty in the linear transfer function estimate is smaller in the region of the spectrum where the energy from the JONSWAP spectrum is located.
The next step of this study will be the application of the methodology to estimate the quadratic transfer function of the same INO WINDMOOR 12MW floater.
below shows a render and scale model of the floater.

Figure 2 .
Figure 2. A JONSWAP spectrum used to generate the data (left).A wave elevation profile sample with the corresponding 1st order force exerted on the floater (right)

Figure 3 .
Figure 3.Comparison of the NARX predictions against the reference value (left).Magnified region from the left graphic (right)

Figure 4 .
Figure 4. Amplitude (left) and phase (right) of the transfer functions obtained by harmonic probing vs. the Hydrostar reference the output of the model for the current time-step n.The lagged values f n−1 , f n−2 , . . ., f n−n f of the load and the lagged values ζ n , ζ n−1 , ζ n−2 , . . ., ζ n−n ζ of the wave elevation profile are concatenated to form the input vector x n