Gaussian mixture autoencoder for uncertainty-aware damage identification in a floating offshore wind turbine

This work proposes an uncertainty-aware approach to the inverse problem of damage identification in a floating offshore wind turbine (FOWT). We design an autoencoder architecture, where the latent space represents the features of the target damage condition. The inverse operator (encoder) is a deep neural network that maps the measurable response to the parameters (means, variances, and weights) of a multivariate Gaussian mixture model. The Gaussian mixture model provides a convenient distributional description that is flexible enough to accommodate complex solution spaces. The decoder receives samples from the Gaussian mixture and maps the damage condition (states) to the system's measurable response. In such a problem, and depending on the quantities being observed (sensor positioning), it is possible that multiple damage states may correspond to similar measurement records. In this context, the main contribution of this work lies in developing a method to quantify the uncertainty within the context of a possibly ill-posed damage identification problem. We employ the Gaussian mixture to express the multimodal solution space and explain the uncertainty in the damage condition estimates. We design and validate the methodology using synthetic data from a FOWT in the commonly adopted OpenFAST software and consider two damage types frequently occurring in mooring lines: biofouling and anchor displacement. The method allows for the estimation of the damage state while capturing the uncertainty in the estimations and the multimodality of the solution under the availability of a limited number of response measurements.

1 Introduction

Floating offshore wind (FOW) is rapidly emerging as a leading form of renewable energy. Faster and steadier winds and a much larger area for future deployment are only a few appealing aspects that make FOW intriguing for industry experts and academics. While FOW's current global capacity barely exceeds 180 MW, recent reports predict that FOW will generate over 250 GW by 2050 (GWEC, 2023). One of the critical barriers limiting the commercial viability of floating offshore wind turbines (FOWTs) is their operational expense, which relates to inspection, monitoring, and faulty component replacement. These tasks are troublesome and expensive in offshore environments, whose access is compromised by weather conditions (Nava et al., 2022).

One of the most critical components of FOWTs is the mooring system. Mooring systems, which inherit the characteristics of similar devices conceived for the oil and gas industry (Årdal et al., 2014), employ steel chains, wire, and cable spanning dozens of meters to anchor the platforms to the seabed or to connect them. They operate under a broad range of environmental and operational conditions that may lead to damage due to corrosion, wear, and fatigue (Li et al., 2018; Liu et al., 2020). These phenomena affect the system's behavior and may compromise the integrity of the platform if no action is taken. Ensuring the safe and optimized performance of mooring lines is crucial to minimizing operational costs and maximizing the profitability of FOWTs.

Condition monitoring of mooring systems is most often carried out on-site through visual underwater inspections (Martinez-Luengo et al., 2016) and ultrasonic testing (Thibbotuwa et al., 2022). Both methods are costly and often inefficient, owing to the requirement for engaging a qualified workforce and appropriate equipment. Companies must also wait for acceptable weather conditions before deploying their crews on FOWT farms to comply with safety requirements. For these reasons, structural health monitoring (SHM) techniques have emerged as a solution to provide a continuous, efficient, and remote assessment of these assets (Ciuriuc et al., 2022; Liu et al., 2023). SHM typically targets an inverse problem solution that aims to identify the condition of a target system and possibly characterize associated damage based on the indirect information delivered via measurements from an instrumentation system (Farrar and Worden, 2013). Two main approaches exist in the field of SHM for FOWTs: model-based (or physics-based, also known as hybrid) and data-based schemes (Liang et al., 2024; Liu et al., 2022b). On the one hand, physics-based models of dynamic systems often employ complex ordinary or partial differential equations (PDEs) that govern the physical phenomena under study. While they require a deeper insight into the underlying physics, these techniques can achieve higher accuracy and generalization at the expense of computational effort (Jonkman, 2007; Wang, 2015; Hall and Goupee, 2015).

On the other hand, data-driven techniques employ extensive datasets to fit the desired outcome (Martinez-Luengo et al., 2016). They can infer highly complex, nonlinear relations, provided that these are witnessed in the available data adopted for training purposes. In the context of SHM for mooring systems, indirect response sensors such as gyroscopes, inclinometers, and GPS trackers have become increasingly attractive (Gorostidi et al., 2023; Coraddu et al., 2024). These devices are relatively inexpensive and easy to deploy while being sensitive to the presence of damage. Experimental data are noisy, constrained to practically measurable quantities, and often limited to a specific condition (e.g., the healthy state), which covers only a subset of the inverse problem's solution subspace. Simulations from a computational parametrization (e.g., a finite-element model) are often employed to complement experimental data and overcome this scarcity (Figueiredo et al., 2010; Zhang and Sun, 2021; Fernandez-Navamuel et al., 2023).

In the past few decades, machine learning (ML) algorithms have gained popularity thanks to advances in data acquisition and transmission, informatics, and computational resources. Particularly interesting are deep neural networks (DNNs), which present advantageous properties such as satisfying the theorem of universal approximation (Hornik et al., 1989) and enabling the incorporation of physical knowledge (Rojas et al., 2024). One key challenge in the resolution of inverse problems is handling uncertainties. Many authors have already dealt with uncertainty when solving inverse problems in other fields. In civil infrastructure, Betancourt et al. (2021) implemented a deep interval neural network to classify damage in a benchmark bridge. Huang and Beck (2015) tested the performance of a sparse Bayesian probabilistic model against incomplete data. Teimouri et al. (2017) designed a Gaussian process-based method to monitor the integrity of a composite airfoil structure. In geophysics, Alyaev and Elsheikh (2022) implemented a mixture density network for the inversion of gamma-ray logs. Liu et al. (2022a) estimated subsurface rock and fluid properties using deep variational autoencoders. Rodriguez et al. (2023) extended the loss formulation stated in Shahriari et al. (2021) and implemented a multimodal variational autoencoder (MVAE) to identify subsurface material properties. In the field of wind energy, Mylonas et al. (2021) proposed a conditional variational autoencoder (CVAE) to deliver uncertainty-robust long-term fatigue predictions in a wind turbine blade based on supervisory, control, and data acquisition (SCADA) signals. Mclean et al. (2023) employed Gaussian processes (GPs) to account for the uncertainty present in power curve models as the damage-sensitive features. Other works propose the use of hierarchical sparse Bayesian learning to solve the system identification problem through model updating using Gibbs sampling (Huang et al., 2017b, a). The authors validated their strategy numerically and experimentally using a relatively simple structure with limited degrees of freedom (DOFs).

Focusing on the application of SHM to assess FOWTs, a vast amount of literature exists using DNNs. For example, Chung et al. (2020) fed response amplitude operator (RAO) data into a DNN to detect anomalies in the cross-section of mooring lines for a tension leg platform. Their research was continued by Lee et al. (2021), who extended their approach to catenary lines and taut mooring systems. Janas et al. (2021) developed a condition-agnostic convolutional neural network (CNN) to detect anomalies caused by the loss of one line in a floating oil and gas vessel, training their model with images of its horizontal displacement history. The recent work by Sharma and Nava (2024) combined CNNs and auto-regressive (AR) models to detect biological fouling-, corrosion-, and anchor-shift-related damage in a floating platform's mooring lines. Their study obtained AR coefficient matrices from displacements and rotations for the DeepCWind OC4 platform's surge, heave, and pitch responses fed as image inputs to the CNN.

Still, many existing works employ deterministic approaches which suffer from significant limitations when tackling real-life inverse problems, such as damage identification via monitoring data. Long-term instrumentation systems are often cheap, very limited, and provide incomplete and noisy measurements. This context makes the solution to the inverse problem highly non-unique and unstable (extremely sensitive to slight changes in the input data) (Adler and Öktem, 2017). Yet in the deterministic scope, Shahriari et al. (2021) proposed a way to define the loss function of DNNs in the measurement space rather than in the ill-posed solution space via an encoder–decoder architecture and a two-step training phase. With this strategy, they constrained the training and prevented undesired solutions. Based on this idea, Gorostidi et al. (2023) attempted to detect failures in the mooring system of a FOWT based on response measurements. They employed statistics-based features from 6 DOFs to identify the level of biofouling and anchoring damage using synthetic data from a FOWT simulated in OpenFAST (Jonkman et al., 2022). Although this approach aids in identifying a physically plausible scenario, it neglects the multimodality of the solution (i.e., various damage scenarios producing the same measured response). Addressing this ill-posedness requires moving from deterministic to probabilistic approaches to quantify the uncertainty in the solution space and provide more reliable assessments.

However, there is still a considerable gap in implementing uncertainty-aware methods for the condition assessment of mooring systems in FOWTs. This work intends to contribute to this direction. Here, we extend our previous work (Gorostidi et al., 2023, 2022), which aimed to infer the underlying health condition of mooring systems using response measurements, which delivered a deterministic assessment. Following the recent work by Rodriguez et al. (2023), we propose a Bayesian approach for quantifying the uncertainty in the delivered damage estimates. Our proposal adapts the MVAE methodology introduced in Goh et al. (2021) and Rodriguez et al. (2023) to solve the inverse problem of damage identification of FOWT mooring systems. The core of the proposed methodology is to probabilistically describe the solution to an inverse problem that may exhibit high multimodality. This situation frequently occurs when dealing with sparse instrumentation systems, which is often the case in SHM applications (Teughels and De Roeck, 2004; Oliveira et al., 2018a, b; Devriendt et al., 2014; Magalhães et al., 2008).

We design an encoder–decoder architecture to address ill-posed inverse problems within the context of damage identification. According to the forward operator, the proposed methodology ensures that the estimated solutions are physically meaningful. We describe the multimodal solution space (damage condition) using a multivariate Gaussian mixture. This parametrized distributional model is mathematically convenient to integrate into the differentiable scheme of DNNs and is sufficiently flexible to accommodate complex distributions. The inverse operator (encoder) estimates the parameters (i.e., means, variances, and weights) that build a corresponding Gaussian mixture describing the distribution of the damage condition features. We draw samples from the estimated posterior distribution model that are subsequently fed into the forward approximation (decoder). With this strategy, we statistically describe the target damage condition space, accounting for the multimodality of the solution, which occurs mainly when the instrumentation system comprises a limited number of sensors.

Training the inverse produces the posterior distribution that describes the damage condition given some input measurements of the system's response. In this work, statistical features of measured rotational DOFs are employed as the response measurements. We employ rotations since these can be measured experimentally through low-cost, low-maintenance sensors, similar to acceleration, while still delivering information that reflects the resulting response. A benefit of using rotations lies in delivering quantities that contain lower-frequency information with respect to acceleration signals. This is useful when the involved response includes lower-frequency components (e.g., rigid-body motions and drifts). The loss function leverages two terms; the first accounts for the measurement misfit, while the second drives the shape of the posterior distribution.

We analyze the effect of uncertainty in the delivered estimates of our proposed scheme in eight test examples that correspond to representative damage scenarios. We further explore how using 1 single DOF increases the ill-posedness of the inverse SHM problem, demonstrating the strength of our proposed method in providing a more reliable diagnosis with poor instrumentation systems. Finally, we further investigate the robustness of the methodology to different measurement error levels. The results demonstrate that our method successfully captures the uncertainty in the predictions, describing the multimodality of the solution mainly in the absence of some response signals.

Despite the successful results, this work suffers certain limitations that ought to be acknowledged. First, the methodology provides a way to describe the uncertainty that is inherent to the ill-posedness of the inverse problem. Such uncertainty induces the multimodality of the output (i.e., the damage condition estimate). The method also reflects the effect of aleatory uncertainty (noisy measurements) as it transfers from the measured data to the estimated condition. However, this work neglects the epistemic uncertainty, which occurs when making a prediction on measurements that correspond to a damage condition that is far from those employed in the training stage. Accounting for such uncertainty and disentangling both sources is beyond the scope of this work and requires further study. Second, since the Gaussian mixture is a parametric approach, it constrains the outcome and prevents a complete characterization of the effect of uncertainty. Improving the distributional model requires many components, which enormously increases the number of parameters to be estimated. This is an important limitation of the method, mostly when scaling up to higher-dimensional spaces. Finally, due to the current scarcity of experimental data from real operating FOWTs, this work is entirely restricted to synthetic data from computational simulations. Integrating experimental data with synthetic scenarios is a key challenge to proving the applicability of the suggested methodology in real-field data.

The remainder of this article is structured as follows. Section 2 derives the mathematical formulation describing the inverse problem. Section 3 presents the turbine–platform assembly and the excitation and damage conditions employed in our simulations. Section 4 describes the specifications of the proposed DNN architecture and training stage. Sections 5 and 6 discuss the method's performance considering (i) 3 DOFs and (ii) a single DOF. Finally, Sect. 7 highlights the conclusions and limitations of the proposed work and reveals future lines of research.

2 Methodology

We describe the rigid-body response of a FOWT platform in terms of 6 DOFs: surge (forward–backward motion), sway (sideways motion), heave (vertical motion), roll (rotation about the longitudinal axis), pitch (rotation about the transverse axis), and yaw (rotation about the vertical axis) (Tran and Kim, 2015). Since transmission of time domain signals is extremely expensive, their content is often condensed in the form of statistical features, which is the typical approach to storing SCADA data (Gorostidi et al., 2023). Let m∈ℳ denote the platform's response with M features extracted from the time domain signals of the 6 DOFs. This response relates to the loading and the system properties through a set of partial differential equations (PDEs). These PDEs describe the aerodynamics, hydrodynamics, servodynamics, and elastodynamics of the coupled system operating under wind and wave excitation conditions w∈𝒲 (Tran and Kim, 2015).

Damage in the mooring system affects its physical properties (e.g., stiffness). We denote the set of damage features describing the condition of the FOWT mooring system with z∈𝒵. These features indicate the level of each possible existing damage and succinctly represent the changes in the coefficients of the governing PDEs. We define $\mathcal{F}:\mathcal{Z}\times \mathcal{W}\to \mathcal{M}$ as the forward operator mapping the damage and loading conditions to the response of the system. Here, ℱ includes the PDEs governing the response of the floating platform of the FOWT. The Cartesian product of the domains 𝒵 and 𝒲, represents all combinations of damage and loading conditions.

2.1 Deterministic inverse operator

In damage identification, we seek the inverse operator $\mathcal{I}:\mathcal{M}\times \mathcal{W}\to \mathcal{Z}$, which, given some noisy response measurements m under a prescribed excitation w, yields the system's damage condition z. This relationship is often unknown and highly nonlinear. Here, we approximate ℐ using deep neural networks (DNNs). DNNs have remarkable power in approximating complex nonlinear functions (Hornik, 1991). Moreover, once trained, such models produce evaluations in milliseconds. These particularities render DNNs appealing in many areas, including damage identification. Let ℐ_θ be the DNN described by the parameters θ that approximates the inverse problem ℐ. For a certain observed input [m,w], we define a loss function to measure the discrepancy between the estimates of the approximate inverse ℐ_θ and the true damage condition of the system, z:

$$\begin{array}{}\text{(1)}& {\mathcal{L}}_{\mathcal{Z}}\left(\mathit{\theta }\right)=\Vert \mathbf{z}-{\mathcal{I}}_{\mathit{\theta }}\left(\right[\mathbf{m},\mathbf{w}\left]\right){\Vert }_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

where we employ the squared l₂ norm as the discrepancy metric. We can find the optimal parameter set θ by minimizing the loss function in Eq. (1) over a training dataset that contains D labeled responses $\mathcal{D}=\mathit{\{}{\mathbf{m}}_i,{\mathbf{w}}_i,{\mathbf{z}}_i{\mathit{\}}}_{i=\mathrm{1}}^D$.

However, in real physical problems, the instrumentation system is often sparse and thus not fully sensitive to the targeted damage scenarios. This situation renders the forward operator ℱ non-injective; i.e., more than one system state (damage condition) may produce the same response (Tarantola, 2004). Consequently, the inverse problem is ill-posed, meaning it might have multiple solutions for the same noisy input. Given this ill-posedness of the inverse problem under incomplete and noisy data availability, minimizing ℒ_𝒵(θ) may produce an infeasible outcome that simply results by averaging possible candidate solutions. Thus, defining the loss on the space of the damage conditions 𝒵 is inconvenient. To overcome this issue, we adopt an encoder–decoder strategy, where the encoder approximates the inverse problem (given the measurements m, estimate the damage condition features z), and the decoder corresponds to the forward operator. With this architecture, the damage condition estimates are enforced to satisfy the forward operator (governing ordinary or partial differential equations) and thus are consistent with the underlying physics describing the system's behavior. This composition of the forward with the inverse, which constitutes the identity mapping, enables the loss function to be expressed on the space of the measured responses, ℳ:

$$\begin{array}{}\text{(2)}& {\mathcal{L}}_{\mathcal{M}}\left(\mathit{\theta }\right)=\Vert \mathbf{m}-\mathcal{F}\circ {\mathcal{I}}_{\mathit{\theta }}\left(\right[\mathbf{m},\mathbf{w}\left]\right){\Vert }_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

using the l₂ norm. Minimization of this loss function ensures that the DNN will report one out of all feasible solutions.

The main bottleneck when minimizing the loss function in Eq. (2) owes to the effort required to massively evaluate the forward operator ℱ. We define a DNN ℱ_ϕ described by the parameters ϕ that approximates the forward operator ℱ. Once trained, ℱ_ϕ provides computationally efficient system responses that substitute the expensive forward evaluations. Figure 1 graphically describes the connection of the forward with the inverse operator using a fully connected architecture.

Figure 1 DNN architecture connecting the forward ℱ_ϕ with the inverse ℐ_θ. The output of the entire DNN yields the reconstruction of the input measurements, ${\widehat{\mathbf{m}}}_i={\mathcal{F}}_{\mathit{\varphi }}\circ {\mathcal{I}}_{\mathit{\theta }}\left({\mathbf{m}}_i\right)$. This architecture permits the definition of the loss function in the measurements space, ℳ. The wind and wave excitations w are fed to both the forward and the inverse. For representation feasibility, we depict the distribution as a bi-variate Gaussian mixture.

We adopt a two-step training strategy proposed and employed in previous works (Shahriari et al., 2021; Rodriguez et al., 2023; Gorostidi et al., 2023). The contribution of these works and their connection with the present approach are extensively described in Sect. 1. We first obtain the optimal forward parameters ϕ^∗:

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}^{\ast }:=\underset{\mathit{\varphi }}{\text{arg\hspace{0.17em}min}}\Vert \mathbf{m}-{\mathcal{F}}_{\mathit{\varphi }}\left(\right[\mathbf{z},\mathbf{w}\left]\right){\Vert }_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

Then, we find the optimal inverse parameters θ^∗:

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}^{\ast }:=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}min}}\Vert \mathbf{m}-({\mathcal{F}}_{\mathit{\varphi }}^{\ast }\circ {\mathcal{I}}_{\mathit{\theta }})\left(\right[\mathbf{m},\mathbf{w}\left]\right){\Vert }_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

The second step incorporates ${\mathcal{F}}_{\mathit{\varphi }}^{\ast }$ as a non-trainable architecture.

After the two-step training, we essentially are interested in using ${\mathcal{I}}_{{\mathit{\theta }}^{\ast }}$ to estimate the damage condition z of the FOWT from the measured responses m and operational conditions w. In this deterministic framework, the inverse operator outputs a single-value estimate for each damage condition given an input [m,w]. This approach precludes adequate interpretation of the uncertainty in the estimates, which are assumed to be 100 % confident. However, despite ensuring that the provided estimate corresponds to one of the possible solutions (assuming an adequate and successful training process), this might differ from the actual condition of the system.

2.2 Bayesian inverse operator

In this work, we overcome the limitation of the deterministic solver by adopting a Bayesian approach to approximate the inverse operator ℐ (Rodriguez et al., 2023; Goh et al., 2021). Since the physical law governing the system is deterministic, we consider the forward to be a deterministic operator. We assume that uncertainty resulting from recording noisy measurements propagates from the data to the predicted damage condition estimates (through the inverse operator). Hence, instead of mapping the measurements to the damage condition, we seek the operator that maps the input probability distribution (of the measurements) to the distribution of the damage condition. The distribution of the damage state features is conditioned on the measured data; i.e., we search for a conditional probability distribution.

Let us assume that the measured response, loading conditions, and unknown system's damage conditions are represented as random variables. Vectors m, w, and z correspond to realizations of each of these sets of random variables. Considering the existence of additive noise that is inherent to the data acquisition process, we express the forward problem as

$$\begin{array}{}\text{(5)}& \mathbf{m}=\mathcal{F}(\mathbf{z},\mathbf{w})+\mathit{ϵ},\end{array}$$

where we substitute ℱ by its optimal approximation ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$, and ϵ is the unknown measurement error with known statistics described by a probability density function (PDF) ϵ∼p(ϵ). Although z is unknown, we can represent its uncertainty using a conditional probability distribution considering its relationship with the measured variables, $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. For any damage condition z, the conditional PDF $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ is the target posterior distribution that forms the estimated target of the inverse problem ℐ under the use of a Bayesian approach. We employ Bayes' theorem to express it with the following proportionality:

$$\begin{array}{}\text{(6)}& p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})\propto p(\mathbf{m}|\mathbf{z},\mathbf{w})\cdot p\left(\mathbf{z}\right),\end{array}$$

where p(z) represents the prior PDF of the unknown damage condition z, and $p\left(\mathbf{m}\right|\mathbf{z},\mathbf{w})$ is the likelihood model that expresses the interrelation between the measurements and the damage condition.

Due to the intractability of the true posterior (usually an unknown non-parametric PDF), we define an inverse operator ℐ_θ that estimates the parameters of an approximate PDF for each damage condition in z, given the response measurements m and operating loads w. Since we have no prior knowledge, we approximate this PDF using a flexible one, given by a mixture of multivariate Gaussian functions (Deisenroth et al., 2020), such that

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{q}_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})\sim \text{GMM}(\mathit{\mu },\mathit{\sigma })={\sum }_{k=\mathrm{1}}^K{\mathit{\alpha }}_k\mathcal{N}(\mathbf{z}|{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k),\\ {\displaystyle }& {\displaystyle }\begin{array}{rl}\mathcal{N}\left(\mathbf{z}\right|{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k)=& {\displaystyle \frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{)}^{Z/\mathrm{2}}\right|{\mathrm{\Sigma }}_k{|}^{\mathrm{1}/\mathrm{2}}}}\\ & \times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{z}-{\mathit{\mu }}_k){\mathrm{\Sigma }}_k^{-\mathrm{1}}(\mathbf{z}-{\mathit{\mu }}_k)\right),\end{array}\end{array}\end{array}$$

where the PDF $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ is a Gaussian mixture model (GMM) dependent on the parameters θ of the DNN ℐ_θ, α_k is the weight of the kth Gaussian in the mixture, and $\mathcal{N}\left(\mathbf{z}\right|{\mathit{\mu }}_k,{\mathrm{\Sigma }}_k)$ is the corresponding multivariate Gaussian distribution with mean vector μ_k and diagonal covariance matrix Σ_k=diag(σ_k) for $k=\mathrm{1},\mathrm{\dots },K$ Gaussians. Compared to the deterministic approach, instead of a single value estimate, now the output of the inverse operator ℐ_θ produces the set of properties describing the mixture, namely the vector ${\text{GMM}}_{\text{props}}=\mathit{\{}{\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\mu }}_K,{\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_K,{\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_K\mathit{\}}$. This results in an output of $K(\mathrm{2}\times Z+\mathrm{1})$ dimensions. For simplicity in notation, we have omitted the dependency of the GMM properties on the DNN parameters θ.

Realizations from $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ represent samples of the damage condition z that are likely produced by the unknown true posterior $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. All the realizations $\mathit{\{}{\mathbf{z}}_h{\mathit{\}}}_{h=\mathrm{1}}^H$ share the same operating conditions w. Recovering the strategy of composing the forward with the inverse (see Fig. 1), we feed the samples to the optimal forward operator ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$ trained in Eq. (3). Figure 2 schematically represents the architecture in the Bayesian approach.

Figure 2 DNN architecture in the adopted Bayesian approach. The output of the inverse operator produces the parameters that describe the PDF of the damage condition z. A sampling layer draws H random samples from the distribution, which are then fed to the optimal forward operator ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$. The output of the entire DNN yields the reconstruction of the input measurements, ${\widehat{\mathbf{m}}}_i={\mathcal{F}}_{\mathit{\varphi }}\circ {\mathcal{I}}_{\mathit{\theta }}\left({\mathbf{m}}_i\right)$.

We aim to minimize the discrepancy between the true posterior $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ and the approximate posterior $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ obtained from ℐ_θ. We first assume that noise ϵ is mutually independent with respect to the unknown damage condition z (Goh et al., 2021). Thus, according to Eq. (5), we can express the likelihood as

$$\begin{array}{}\text{(8)}& p\left(\mathbf{m}\right|\mathbf{z},\mathbf{w})=p(\mathbf{m}-\mathcal{F}\left(\right[\mathbf{z},\mathbf{w}\left]\right)).\end{array}$$

We then assume the noise follows a Gaussian distribution, $p\left(\mathit{ϵ}\right)=\mathcal{N}(\mathrm{0},\mathrm{\Gamma })$, where $\mathrm{\Gamma }=\text{diag}\left(\mathit{\beta }\mathcal{F}\right([\mathbf{z},\mathbf{w}]){)}^{\mathrm{2}}$ is a vector that contains the non-zero elements of a diagonal matrix, and the parameter β corresponds to the noise level. This allows us to rewrite Eq. (6) as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}p& \left(\mathbf{z}\right|\mathbf{m},\mathbf{w})\propto p(\mathbf{m}-\mathcal{F}(\mathbf{z},\mathbf{w}))\cdot p(\mathbf{z})={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{)}^{M/\mathrm{2}}\right|\mathrm{\Gamma }{|}^{\mathrm{1}/\mathrm{2}}}}\\ & \times \exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{m}-\mathcal{F}([\mathbf{z},\mathbf{w}]\left){)}^t{\mathrm{\Gamma }}^{-\mathrm{1}}\right(\mathbf{m}-\mathcal{F}\left(\right[\mathbf{z},\mathbf{w}\left]\right))\right),\end{array}\end{array}$$

where p(z) is the prior distribution of the damage condition properties, which can follow any PDF. The Kullback–Leibler divergence (KL) is defined as a statistical measure of the distance between two PDFs (Deisenroth et al., 2020; Asperti and Trentin, 2020). KL assumes a null value when the compared PDFs are equal and can be defined as

$$\begin{array}{}\text{(10)}& \text{KL}\left[p\right(x\left)\right|\left|q\right(x\left)\right]=\int p\left(x\right)\log {\displaystyle \frac{p\left(x\right)}{q\left(x\right)}}\mathrm{d}x,\end{array}$$

where x indicates a realization of the random variable X, and p(x) and q(x) are the two different PDFs. Here, we use the KL metric to evaluate the distance between the true posterior $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ and its approximation estimated by the inverse operator ℐ_θ, denoted as $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$, yielding

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\text{KL}& \left[q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})\left|\right|p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\right]\\ & =\int q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})\log {\displaystyle \frac{q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})}{p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})}}\mathrm{d}z.\end{array}\end{array}$$

However, the KL divergence term exhibits certain shortcomings that weaken its strength as a distance metric; it is asymmetric, it does not satisfy the triangle inequality, and it produces an intractable term (the evidence of the data distribution p(m)) (Blei et al., 2017). Instead, a lower bound is calculated for the evidence, known as the evidence lower bound (ELBO) (Blei et al., 2017). ELBO is the loss function commonly employed in variational autoencoders (VAEs) to account for the discrepancy between the two distributions (Goh et al., 2021; Rodriguez et al., 2023). We obtain the ELBO loss by exploiting the KL expression and isolating the intractable terms:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\text{KL}[q_{\mathit{\theta }}& \left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\right|\left|p\right(\mathbf{z}|\mathbf{m},\mathbf{w})]\\ =& \int q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})\log {\displaystyle \frac{q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})}{p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})}}\mathrm{d}z\\ =& \int q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})\mathrm{d}z\\ & -\int q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\log p\right(\mathbf{z}|\mathbf{m},\mathbf{w})\mathrm{d}z\\ =& {\mathbb{E}}_q\left[\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})]-{\mathbb{E}}_q[\log p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\right]\\ =& {\mathbb{E}}_q\left[\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})]-{\mathbb{E}}_q\left[\log {\displaystyle \frac{p(\mathbf{m},\mathbf{w},\mathbf{z})}{p(\mathbf{m},\mathbf{w})}}\right]\\ =& {\mathbb{E}}_q\left[\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})-\log p(\mathbf{m},\mathbf{w},\mathbf{z})]\\ & +\log p(\mathbf{m},\mathbf{w}).\end{array}\end{array}$$

In Eq. (12), we employ the relationship between the joint distribution and the posterior: $p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})=\frac{p(\mathbf{m},\mathbf{w},\mathbf{z})}{p(\mathbf{m},\mathbf{w})}$. We rearrange the terms and define the ELBO as

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\text{ELBO}=& \log p(\mathbf{m},\mathbf{w})-\text{KL}\left[q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})\left|\right|p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\right]\\ =& {\mathbb{E}}_q\left[\log p\right(\mathbf{m},\mathbf{w},\mathbf{z})-\log q_{\mathit{\theta }}(\mathbf{z}|\mathbf{m},\mathbf{w})]\\ =& {\mathbb{E}}_q\left[\log \right(p\left(\mathbf{m}\right|\mathbf{w},\mathbf{z}\left)p\right(\mathbf{w}\left)p\right(\mathbf{z}\left)\right)]\\ & -{\mathbb{E}}_q\left[\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})]\\ =& {\mathbb{E}}_q\left[\log p\right(\mathbf{m},\mathbf{w}\left|\mathbf{z}\right)]+{\mathbb{E}}_q[\log p\left(\mathbf{z}\right)]\\ & -{\mathbb{E}}_q\left[\log q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})],\end{array}\end{array}$$

where we assume that the operating conditions w and the damage properties z are independent (i.e., $p(\mathbf{w},\mathbf{z})=p\left(\mathbf{w}\right)\cdot p\left(\mathbf{z}\right)$). Hence, we can remove p(w) in the second line of the equation, as it is independent of z.

For a certain observation {m,w}, we draw H samples from the posterior and approximate the ELBO loss function as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathcal{L}}_{\text{ELBO}}\left(\mathit{\theta }\right)\approx \\ & {\displaystyle \frac{\mathrm{1}}{H}}\sum _{h=\mathrm{1}}^H[\underset{\text{Likelihood}}{\underbrace{\log p(\mathbf{m},\mathbf{w}|{\mathbf{z}}^h)}}+\underset{\text{Prior}}{\underbrace{\log p\left({\mathbf{z}}^{\mathbf{h}}\right)}}-\underset{\text{Approx. posterior}}{\underbrace{\log q_{\mathit{\theta }}\left({\mathbf{z}}^h\right|\mathbf{m},\mathbf{w}\left)\right]}}.\end{array}\end{array}$$

By substituting the likelihood from Eq. (9) into Eq. (14), we finally express ℒ_ELBO as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\mathcal{L}}_{\text{ELBO}}\left(\mathit{\theta }\right)\approx \\ & \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{H}}\sum _{h=\mathrm{1}}^H[& -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{m}-\mathcal{F}([{\mathbf{z}}^h,\mathbf{w}]\left){)}^t{\mathrm{\Gamma }}^{-\mathrm{1}}\right(\mathbf{m}+\mathcal{F}\left(\right[{\mathbf{z}}^h,\mathbf{w}\left]\right))\\ & +\log p\left({\mathbf{z}}^{\mathbf{h}}\right)-\log q_{\mathit{\theta }}\left({\mathbf{z}}^h\right|\mathbf{m},\mathbf{w}\left)\right].\end{array}\end{array}\end{array}$$

Here, the first term accounts for the data misfit, which is the error between the true measurements and the reconstructions provided by ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$. The second term refers to the prior, which we assume to follow a bounded uniform distribution $p\left(\mathbf{z}\right)\sim \mathcal{U}[{\mathbf{b}}_{\text{low}},{\mathbf{b}}_{\text{up}}]$ with lower and upper bounds b_low and b_up, respectively. The last term measures the probability that the hth sample belongs to the estimated distribution $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The second term can be neglected by directly constraining the Gaussian mixture density function to the desired interval (according to the assumption of uniform distribution). Note that minimizing the Kullback–Leibler divergence between the original and the estimated posteriors ($\text{KL}\left[q_{\mathit{\theta }}\right(\mathbf{z}|\mathbf{m},\mathbf{w})\left|\right|p\left(\mathbf{z}\right|\mathbf{m},\mathbf{w}\left)\right]$) is equivalent to maximizing ℒ_ELBO or minimizing its negative. For a training dataset 𝒟 with D labeled observations ($\mathcal{D}=\mathit{\{}{\mathbf{m}}_i,{\mathbf{w}}_i,{\mathbf{z}}_i{\mathit{\}}}_{i=\mathrm{1}}^D$) and H drawn samples from the estimated PDF, we obtain the optimal parameter set θ^∗ by minimizing the negative of ℒ_ELBO:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathit{\theta }}^{\ast }& :=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}min}}{\displaystyle \frac{\mathrm{1}}{D\cdot H}}\\ & \begin{array}{rl}\times \sum _{i=\mathrm{1}}^D\sum _{h=\mathrm{1}}^H[& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathbf{m}}_i-\mathcal{F}([{\mathbf{z}}_i^h,{\mathbf{w}}_i]\left){)}^t{\mathrm{\Gamma }}^{-\mathrm{1}}\right({\mathbf{m}}_i-\mathcal{F}\left(\right[{\mathbf{z}}_i^h,{\mathbf{w}}_i\left]\right))\\ & +\log q_{\mathit{\theta }}\left({\mathbf{z}}_i^h\right|{\mathbf{m}}_i,{\mathbf{w}}_i\left)\right].\end{array}\end{array}\end{array}$$

3 Case study

We consider the use case of the 5 MW FOWT designed by the National Renewable Energy Laboratory (NREL) (Jonkman et al., 2009) mounted atop the DeepCWind semisubmersible platform (Robertson et al., 2014a). The model was developed as part of the Offshore Code Comparison Collaboration Continuation (OC4) project, which aimed to verify the accuracy of offshore wind turbine dynamics simulation tools through extensive code-to-code comparisons. This involved contributions from numerous organizations worldwide, ensuring a robust and well-validated design. The semisubmersible floater for OC4 was specifically designed to serve as a benchmark for offshore wind energy research and development purposes (Robertson et al., 2014b). The motivation for adopting such a turbine–floater assembly in this study lies in its computational availability and popularity within the FOW research community. Moreover, the OC4 model underwent experimental validation. The validation process included comparisons between numerical simulations and experimental data from wave tank tests (Borisade et al., 2018). The OC4 floater, depicted in Fig. 3, provides stability to the turbine facing unsteady and unpredictable sea conditions thanks to three partially ballasted base cylindrical columns. A set of pontoons and cross-braces connects each column to the others and to the main tower. The platform is held in place by three catenary mooring lines located 120° from one another. The selected measurements describing the platform's response include features from the 3 rotational DOFs: roll, pitch, and yaw. From a sensory perspective, these measurements can be experimentally obtained using inclinometers, cheap and reliable devices often used in long-term monitoring, or via micro-electromechanical system (MEMS) technologies, given the mechanical symmetry of the system in terms of force and load distribution. Compared to acceleration signals, which are often measured in the field of FOWT condition assessment, displacements (and coupled-rotation DOFs) tend to be more sensitive under low-frequency dynamics, such as floating-platform responses. In the recent work by Sharma and Nava (2024), the authors analyzed the sensitivity of acceleration and displacement signals under mooring system damage (anchoring), highlighting the potential of displacements over accelerations even for low-level damage.

Figure 3 5 MW NREL FOWT mounted atop the DeepCWind OC4 semisubmersible platform.

Figure 4 describes the response under common environmental conditions. We extract a set of modal statistics from these responses, which we assume are sufficiently descriptive of the platform's movement.

Figure 4 Response amplitude operator (RAO) of the FOWT for roll response in (a) time and (b) frequency domain, for pitch response in (c) time and (d) frequency domain, and for yaw response in (e) time and (f) frequency domain.

As observed in Fig. 4, owing to the system's symmetry, rotations in the roll DOF are harmonic oscillations around the mean position with an energy content associated with the natural surge and roll periods, while pitch exhibits a more colorful power spectral density (PSD) being affected by the external excitation.

From the time domain responses (see Fig. 4a, c, and e), we compute the mean displacement as

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{‾}}{x}=\underset{t_{\mathrm{0}}}{\overset{t_f}{\int }}x\mathrm{d}t\approx {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{x}_i,\end{array}$$

where N indicates the total number of data points. We further compute the standard deviation of the response as

$$\begin{array}{}\text{(18)}& \mathit{\sigma }=\sqrt{{\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{i=\mathrm{1}}^N(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}.\end{array}$$

We assume the time domain response to be stationary by neglecting its transient state. To obtain the frequency spectra (see Fig. 4b, d, and f), we compute the PSD of the signals, which describes the power distribution across a frequency range as (Theodoridis, 2020)

$$\begin{array}{}\text{(19)}& S_x\left(f\right)=\underset{T\to \mathrm{\infty }}{lim}{\displaystyle \frac{\mathbb{E}\left[\right|F_x\left(f\right){|}^{\mathrm{2}}]}{\mathrm{2}T}},\end{array}$$

where F_x(f) is the Fourier transform of the time domain signal for any DOF x. Finally, we identify two dominant peak frequencies as

$$\begin{array}{}\text{(20)}& {\displaystyle }{f}_{\mathrm{1}}=\underset{f\in \left[\mathrm{0},f_{\text{thresh}}\right]}{\text{arg\hspace{0.17em}max}}{S}_x\left(f\right),\text{(21)}& {\displaystyle }{f}_{\mathrm{2}}=\underset{f\in \left[f_{\text{thresh}},\mathrm{\infty }\right]}{\text{arg\hspace{0.17em}max}}{S}_x\left(f\right).\end{array}$$

One of the peaks for each DOF usually matches the platform's natural frequency. In contrast, the further spectral peaks reflect the influence of external conditions, e.g., wind, current, and waves, in the system's response (Benitz et al., 2014). We have split each DOF's natural and excitation frequencies using threshold frequencies f_thresh to ensure both peaks are identified.

To assess the magnitude of the peaks and the intensity of all the frequencies in the spectra, we also measure the zero-th momentum as (Sundar, 2017)

$$\begin{array}{}\text{(22)}& m_{\mathrm{0}}=\underset{\mathrm{0}}{\overset{{\mathit{\omega }}_f}{\int }}\mathit{\omega }{S}_x\mathrm{d}\mathit{\omega },\end{array}$$

where ω is the angular frequency in radians per second. These are the five features we employ to describe the time domain response of the platform. These statistics are ultimately the inputs of our neural network. This method allows for a significant reduction in the number of inputs required to assess the health status of mooring lines' integrity while maintaining high physical accuracy.

In this work, we use NREL's open-source wind turbine simulation tool OpenFAST (OpenFAST Documentation, 2023), which evaluates the influence of aerodynamic (Jonkman et al., 2015) and hydrodynamic (Jonkman et al., 2014) excitations on the response of the floating platform. The number of deployed real-scale FOW turbines is limited; moreover, the data related to the performances of these devices are proprietary, and even if they were publicly available, the amount of labeled data under damage conditions in the mooring systems might be extremely reduced, given the short life of these platforms. The current lack of operational data from such devices demands solutions that may exploit simulation tools in the delivery of predictive models. OpenFAST is regarded as a highly accurate and reliable tool for numerical simulations of FOWTs in wind–wave environments (Yang et al., 2021; Reig et al., 2024; Rinker et al., 2020) owing to existing validation efforts using experimental data (Coulling et al., 2013; Robertson et al., 2017). The OC4 model, in particular, has been validated in OpenFAST against experimental data for free-decay tests (Gorostidi et al., 2023; Liu et al., 2019).

Our simulations aim at solving the system of equations given by

$$\begin{array}{}\text{(23)}& (M+A_{\mathrm{\infty }})\ddot{q}+Kq=\sum F(t,\mathit{\omega }),\end{array}$$

where M is the system's mass matrix, A_∞ is the added mass matrix, and K is the hydrostatic stiffness matrix. Terms q and $\ddot{q}$ encompass the system's position and acceleration, respectively (Faltinsen, 1993). The force term can be split into components reflecting the contribution of wave, wind, viscous, mooring, and radiation-damping forces (Jonkman and Matha, 2011).

The simulation process follows the same structure as that presented by Gorostidi et al. (2023). We assign the environmental conditions for each simulation by selecting a combination for significant wave height $H_{\mathrm{S}}\in [\mathrm{2},\mathrm{15}]$ (m) and peak period $T_{\mathrm{P}}\in [\mathrm{1},\mathrm{15}]$ (s). The decision to use these features is motivated by their common availability in real practice since supervisory control and data acquisition (SCADA) systems typically include such measurements.

We have defined evenly spaced values for both H_S and T_P within their feasible interval, and one combination is randomly selected for each simulation using Monte Carlo sampling. These two variables define a Pierson–Moskowitz spectrum, which estimates the distribution of the energy of ocean waves based on their frequency using the empirical correlation proposed by Pierson and Moskowitz (1964). The integration of this spectrum defines the temporal evolution of the wave force component of the total force in Eq. (23). We select wind velocity W_V in a similar manner, with speeds ranging from 1 to 30 m s⁻¹. In this work, we consider uniform wind speed profiles.

We then introduce damage to the mooring lines in MoorDyn, OpenFAST's mooring line dynamics module (Hall, 2020). In this work, as shown in Fig. 5, we distort one of the lateral mooring lines of the platform.

Figure 5 Top view of the DeepCWind OC4 semisubmersible platform. The mooring line exposed to damage is highlighted in red.

Most failures in the mooring system occur during the operational phase, according to a survey focused on floating production storage and offloading (FPSO) platforms for the oil and gas sector (Fontaine et al., 2014). Pitting corrosion, fatigue due to cyclic loading, and abrasion with the seabed may represent some of the most frequent causes leading to critical failure of the mooring systems. In the present study, we analyze the effect of degradation caused by two other common forms of damage: biological fouling (Decurey et al., 2020) and anchor point slippage (Liu et al., 2021; Sharma and Nava, 2024). These damage mechanisms affect the mechanical properties of the platform's mooring lines, e.g., mass, stiffness, and buoyancy, and may accelerate wear or cause premature failure (Spraul et al., 2017). Biofouling, in particular, is a slow process that affects the mass and drag of the mooring system, affecting the performance and stability of the platform. Anchor slipping, moreover, drastically affects the stiffness of the system (Liu et al., 2021).

We simulate biofouling by modifying the mass per unit length and diameter of a segment located at the center of the mooring line. We consider the maximum biofouling damage to increase these properties by 10 %. We induce anchoring damage by displacing the line's anchor points' x and y coordinates, which we assume to have an effect on the line's stiffness, causing alterations in the response of the platform. We consider the maximum anchoring damage at 20 m parallel to its baseline orientation of 240° with respect to the x direction. To label the scenarios, we employ two severity coefficients, one for each damage type, where the maximum damage level corresponds to a value equal to one, and healthy mooring line scenarios are assigned null coefficient values. Any intermediate damage coefficient reflects mild degradation, which modifies the line's properties following a linear interpolation. We sample these coefficients by sampling from a folded Gaussian distribution located around zero, with a standard deviation of 0.35. The rationale behind this is that we seek a training dataset containing mostly low-severity scenarios, much more frequent in common operation (Gorostidi et al., 2023).

For any scenario, we first sample the damage coefficients describing the damage condition. We subsequently use Monte Carlo to sample the loading conditions and assign different sea states. With this approach, we build a large dataset with cases that include milder or more energetic sea states for all the considered damage scenarios. Each scenario involves a simulation recreating 30 min of FOWT dynamics, computing platform responses every 0.025 s, and subsequently extracting the selected features. In summary, any scenario is defined by (i) 3 features describing the loading conditions, namely H_S, T_P, and W_V; (ii) 15 response features, 5 for each of the 3 considered DOFs; and (iii) the 2 damage condition coefficients. We produce 60 000 samples in parallel batches using 120 Intel Xeon (R) E5-2680, 2.70 GHz CPUs (Donostia International Physics Centre, 2022), taking approximately 42 h.

4 Neural network design and training

We employ TensorFlow 2.13 to treat the datasets and train the Bayesian DNN for damage condition assessment (Abadi et al., 2015). We split our dataset into training (𝒟_train), validation (𝒟_val), and testing (𝒟_test), each containing 70 %, 20 %, and 10 % of the total samples, respectively. We then use the MinMax scaler (Scikit Learn, 2024) to constrain the environmental conditions and response features to the interval [0,1] so as to ensure that high-order features do not outweigh lower-order ones in the training process. The rescaling function is based on the training data and applies to the three datasets.

Next, we specify the DNN architecture and the hyperparameters describing the training routine. There are no general guidelines regulating the pursuit of an optimal configuration. Numerous hyperparameters exist in this case, pertaining to the model definition and its training process, including, among others, layer counts, types, and sizes; activation functions; and regularization techniques. This renders the search space virtually infinite. As a common practice, developers usually test the performance of a range of candidate architectures designed based on educated guesses and experience, aiming to strike a balance between computational efficiency and prediction accuracy.

In this work, we employ a combination of hyperbolic tangent (Namin et al., 2009) and rectified linear unit (ReLU) (Agarap, 2018) functions for the hidden layers. We implement the weight initialization method proposed by Aldirany (2024), who suggested that TensorFlow's default Glorot uniform initialization scheme (Glorot and Bengio, 2010) is unsuitable for non-differentiable and non-zero mean functions, such as ReLU. Instead, we apply the He uniform initialization (He et al., 2015) to the ReLU layers. At the output layer, we use three different activations for the properties of the multivariate Gaussian mixture: the sigmoid function (Han and Moraga, 1995) for the means, as we seek a smooth function to estimate damage condition coefficients in the interval [0,1]; the softplus function (Zheng et al., 2015) for the variances, as a smooth equivalent of ReLU to enforce positive values; and the softmax function (Goodfellow et al., 2016) for the weights so that their sum is equal to one, and each value falls within [0,1]. We have observed adequate training performance when using the parameters shown in Table 1. However, other configurations may also provide satisfactory results. For the multivariate Gaussian mixture PDF, we assign k=5 Gaussian components. We have no prior knowledge of the optimal number of Gaussians to reflect the uncertainty in the damage condition properties. The larger the number, the more flexible the mixture will be. However, the number of Gaussians directly affects the number of parameters to be estimated; thus, it hampers the training process. After some trial and error analysis, the authors reached adequate results using five Gaussians. We draw H=10 samples to be fed into ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$. A larger number of samples accelerates convergence at the cost of more time-consuming iterations. Finally, we consider a noise parameter β equal to 0.076, which roughly corresponds to noise levels of up to approximately 8.5 % (Bishop, 2006). The value of β is critical to leverage the contribution of both loss terms in Eq. (16) during training and thus to properly describe the uncertainty in the solutions.

We follow the two-step training procedure described in Sect. 2. We first find the optimal forward operator ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$ by minimizing Eq. (3), as described in our previous work (Gorostidi et al., 2023). We pre-train the decoder to approximate the physical law of the system, since this enables access to the derivative quantities that are needed when training the encoder, according to the loss function in Eq. (16). In this manner, the decoder serves to impose the known physical law as a type of inductive bias. A further benefit of this approach is that by specifying fewer unknowns (only those corresponding to the encoder), we reduce the difficulty of the inference task; in particular, we decrease the number of local minima.

Table 1 Specifications of our Gaussian mixture autoencoder. Tanh: hyperbolic tangent; ReLU: rectified linear unit; GU: Glorot uniform, HU: He uniform.

We present the evolution of the decoder's loss in Fig. 6. The parameters of ${\mathcal{F}}_{{\mathit{\varphi }}^{\ast }}$ are frozen for the next training step. We obtain the optimal inverse operator ${\mathcal{I}}_{{\mathit{\theta }}^{\ast }}$ by minimizing the ELBO loss described in Eq. (16). Figure 7 depicts the evolution of the loss function during training, including the total loss value, and the two participating terms.

Figure 6 Evolution of training and validation losses for the forward operator ${\mathcal{F}}_{{\mathit{\varphi }}^{*}}$.

Figure 7 Evolution of training and validation losses for the 3-DOF case: (a) data misfit and mixture density components of the loss and (b) training and validation ELBO losses.

5 Results with 3 DOFs

This section analyzes the damage identification performance using the test dataset 𝒟_test unseen during training and validation tasks. Here, we select eight damage scenarios from 𝒟_test to visualize the results. The eight data points represent different damage conditions and are summarized in Table 2.

Table 2 Properties of the selected examples.

The deterministic approach employed in Gorostidi et al. (2023) minimized the data misfit (ℒ_ℳ), producing a single value estimate for each input measurement. However, despite being a feasible solution, it might be far from the true one when multiple solutions coexist. For any response measurement m, we can identify the feasible solutions as those producing a reduced value of ℒ_ℳ. Accordingly, these solutions must produce a high probability value in the estimated posterior PDF, $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. Figures 8 and 9 depict the solution space and compare two contour maps within the solution space: (i) the density value over $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ and (ii) the ℒ_ℳ value for the eight test example scenarios. The contour plot of data misfits reveals the existence of multiple damage scenarios producing similar responses (measurements), i.e., producing a small value of ℒ_ℳ. Both contours must ideally be identical for the posterior to show all the feasible solutions. However, the constraints and assumptions imposed for the posterior to be tractable (parametrization of the distribution to a Gaussian mixture, diagonal covariance matrix assumption, etc.) restrict the shape of the PDF. For comparison purposes, the figures include both the ground truth and the deterministic solutions produced in Gorostidi et al. (2023).

Figure 8 First four test examples: the left-hand figures represent the contour plot of the estimated posteriors $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The right-hand figures represent the data misfit (ℒ_ℳ) value.

Figure 9 Last four test examples: the left-hand figures represent the contour plot of the estimated probability density functions $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The right-hand figures represent the data misfit (ℒ_ℳ) value.

In the figures, we have constrained the contour maps of ℒ_ℳ to enable proper visualization of the targeted regions. Despite the shape limitations of the estimated posterior distributions, the resulting contours (left-hand figures) exhibit a clear correspondence with the expected solution space described by ℒ_ℳ. Thus, the contour plot of the estimated posterior PDF $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ describes the uncertainty in the solution space.

We observe that, in most cases, the shape accommodates close to a bi-variate Gaussian distribution, indicating the existence of a dominating mode. These results suggest that the selected sensing system and the extracted features suffice to uniquely identify the damage condition causing the observed measurements. Analyzing the deterministic solution, we observe that it lies close to the true damage condition, indicating acceptable predictions. However, as we explore in Sect. 6, employing an incomplete sensing system (fewer DOFs available) might contribute to the ill-posedness of the inverse and thus result in a much more uncertain damage condition.

In this work, we make two main assumptions in the training dataset to constrain the scope of the analysis. First, the results presented in this section assume the availability of some loading condition features describing the wind and wave excitation. Although SCADA systems are frequently found in FOWTs, such measurements may be inaccessible or be subject to failure. Appendix A describes the results obtained when training the inverse operator without any information regarding the loading conditions.

Second, we particularize the methodology to the damage condition of one single mooring line. As a proof of concept, we assume that only one mooring line may suffer damage at a time. However, in real practice, simultaneous damage may happen in two or more components of the target system (e.g., two mooring lines). Appendix B presents the results obtained for a new case study where we consider anchoring damage co-occurring at two mooring lines.

Given the computational limitations of simulating damage scenarios, a relevant aspect that needs to be considered relates to the capability of the DNN to provide a reliable outcome for a certain scenario that was unseen during training (although of the same nature/type as those used for training). Appendix C analyzes such a situation via two different studies to demonstrate the ability of the distributional model to approximate the true solution and provide information about the uncertainty origin.

5.1 Exploring robustness to noise

One critical issue to tackle when dealing with data from real operative systems is the effect of noise. We represent the measurement error as additive noise with a Gaussian distribution. We assume the covariance matrix to be a diagonal matrix $\mathrm{\Gamma }=\left(\mathit{\beta }\mathcal{F}\right([\mathbf{z},\mathbf{w}]){)}^{\mathrm{2}}$, where the scalar β affects the variances of the noise components. A low value of β indicates almost null noise in the data, and the total loss value is mainly due to the data misfit term. Contrarily, a larger value of β is associated with higher noise levels and an increase in the contribution of the distributional term. We can relate β with the signal-to-noise ratio (SNR), which is given by (see Johnson, 2006)

$$\begin{array}{}\text{(24)}& \text{SNR}=\mathrm{10}{\log }_{\mathrm{10}}{\displaystyle \frac{\mathbb{E}\left(S^{\mathrm{2}}\right)}{\mathbb{E}\left(n^{\mathrm{2}}\right)}},\end{array}$$

where S represents the measured signal and n the noise. For white Gaussian noise with a null mean, we substitute the expectation of the noise by its variance. As an indicative value of the noise level based on the training dataset, assigning β=0.05 corresponds to a SNR≈26 dB.

We explore increasing noise levels by modifying the value of β. Figure 10 depicts this effect for the first test example in Table 2. Analogous results are observed for the other examples. The results reveal the sensitivity of the estimated PDF to the noise level. We observe that for β=0.1 (SNR≈20 dB), the entire solution space is feasible for the biofouling damage represented by z₁. This demonstrates a limited sensitivity of the measured features to this damage compared to anchoring (z₂), which can still be estimated. These results support the robustness of the method in the presence of high noise levels.

Figure 10 Analysis of noise effect in the $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$ results for the first test example. Increasing noise levels according to (a) β=0.025 (SNR≈32 dB), (b) β=0.050 (SNR≈26 dB), (c) β=0.075 (SNR≈22.5 dB), and (d) β=0.10 (SNR≈20 dB).

6 Results with 1 DOF

To highlight the potency of the proposed method, we explore the multimodality of the solution space when the instrumentation system is sparse (limited number of sensors). This situation is very common when low-cost, long-term sensing devices are installed. Measuring accelerations is also an extended and cost-effective practice. Damage in the mooring elements affects the dynamics of the entire FOWT system (including the 6 degrees of freedom). However, the dynamics of the mooring system correspond to low frequencies and strongly exciting sway DOF, which is coupled with the rotation DOF roll. Compared to accelerations, displacements (and coupled-rotation DOFs) tend to be more sensitive to damage in the mooring systems, such as anchor displacements, in terms of signal power (Sharma and Nava, 2024). To explore the multimodality in the damage condition estimates, we have selected the roll DOF as the only available signal since it is particularly sensitive to damage given the mechanical symmetry of the system and the unidirectionally of the external excitation and aids in the visualization of the results.

In this case, we employ only five response features – those associated with roll – as the input data for the DNN. We adapt the architecture of ℐ_θ described in Sect. 4 to accommodate the new input dimensions and repeat the training step to minimize ℒ_ELBO. Figures 11 and 12 show the results for the eight damage scenarios.

Figure 11 First four test examples for 1 DOF: the left-hand figures represent the contour plot of the estimated probability density functions $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The right-hand figures represent the data misfit (ℒ_ℳ) value.

Figure 12 Last four test examples for 1 DOF: the left-hand figures represent the contour plot of the estimated probability density functions $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The right-hand figures represent the data misfit (ℒ_ℳ) value.

The data misfit (ℒ_𝒟) contour maps on the right-hand side of Fig. 12 reveal a much more spread range of plausible damage conditions across the solution space for any given input measurement. The contour plots of the estimated probabilities (left-hand figures) match with the ℒ_𝒟 map. Given that the uncertainty has increased in the solution space (ill-posedness), the results reveal a more clear multimodality in the solution space. Compared to Figs. 8 and 9, where the deterministic solution was pretty close to the ground truth, now we observe a considerable separation. This discrepancy already provides information about the existence of multiple solutions. For example, for the test case 1 (see Fig. 11a and b), we clearly observe how the Gaussians try to adopt the tube-shaped contour map of minimal ℒ_𝒟 values that seem to connect the ground truth and the deterministic solution.

Figures 9–12 qualitatively illustrate the multimodal nature of the solution to our inverse problem. Next, we further quantify the performance of our network using common metrics. Table 3 reports the root mean square error (RMSE) of the estimations using both deterministic and probabilistic models. The table's left half presents the RMSE values for the 3-DOF case, while its right half displays the RMSEs for the case of incomplete instrumentation with the roll DOF only. RMSEs were calculated using 10 samples from the Gaussian mixture autoencoder. Since the deterministic model provides a point prediction, the RMSEs of the deterministic predictions are equivalent to the absolute error between the ground truth and the model's estimation.

Table 3 RMSE values for 3-DOF and 1-DOF cases using both deterministic and Gaussian mixture autoencoder models.

The 3-DOF errors in Table 3 concur with the contour plots presented in Figs. 8 and 9, which proved the selected instrumentation sufficient for the deterministic model to accurately estimate mooring degradation. Table 3 shows that the deterministic model's predictions using 1 DOF are poorer, with our Bayesian approach exhibiting improved behavior, especially in Scenarios 1, 3, and 7. These cases, as shown in Figs. 11 and 12, reflect a highly multimodal solution and prove that the deterministic model does not accurately estimate the true solution but rather captures one of the plausible solutions.

These results support the capacity of the proposed method to provide an uncertainty-aware assessment in the form of the probability value for each pair of damage conditions. In such cases, postprocessing can be applied to extract the relevant information from the multivariate contour map, such as defining extrema or establishing a threshold to select the most likely scenarios. Statistical data compression techniques (e.g., principal component analysis, PCA) could be applied to produce a lower-dimensional representation of the space that enables visualization.

7 Conclusions

This work proposes a Bayesian deep learning strategy to achieve an uncertainty-aware assessment of the condition of the mooring system in a floating offshore wind turbine (FOWT). To remedy this, we employ a mixture of multivariate Gaussians to track how the uncertainty is propagated from the available measurements to the estimated health condition diagnostics, providing more robust and reliable estimates. We test the performance of the method using measurements from 3 degrees of freedom and explore the robustness against increasing noise levels, with successful results. We also analyze the benefits of the method when dealing with sparse sensor scenarios, such as measuring at only 1 degree of freedom, revealing the ability of the method to reveal the multimodal nature of the solution in ill-posed scenarios.

In line with the limitations and challenges stated in Sect. 1, here we determine future research lines and challenges to be addressed. Scaling up the methodology to more complex damage spaces (e.g., simultaneous damage in all mooring lines, different directions of the anchor displacements, corrosion effect) is one of the key challenges to be addressed. Given the intractable dimensionality increase in the number of parameters of the Gaussian mixture, we explore alternate distributional models, such as random fields (Birmpa and Katsoulakis, 2021), Gaussian processes (GPs) (Li et al., 2019), or copulas and vine copulas (Letizia and Tonello, 2025).

We further consider a future research line to focus on more exhaustive feature extraction that includes more refined system identification features, including the transmissibility and transmittance, or autocorrelation features, as they may enhance the damage identification task. Also, to achieve a higher-level assessment, we investigate clustering techniques to separate the different types of uncertainty according to the nature and properties of the existing sources (e.g., multimodality of the estimate, measurement noise, model uncertainties, extrapolation over unknown measurements) (Kamariotis et al., 2024).

Finally, it is a core priority of the authors to validate the proposed methodology for a case study where both experimental and synthetic data can be combined. Experimental data will present important limitations, including a strong concentration around a reduced region of the domain (mostly the undamaged condition and other unlabeled damage scenarios). Hence, it is mandatory to find an adequate strategy to complement these data with synthetic simulations, accounting for the computational limitations in terms of time and resources. In this sense, covering the entire range of possible damage scenarios is computationally prohibitive unless it is carefully tackled. Active learning methods enable an efficient design of the synthetic dataset required for the training stage, as they guide the simulation process to enrich the most poorly characterized regions.

Finally, given the potential of the proposed methodology to be applied to different inverse problems, the authors are currently exploring its application to the field of damage identification for bridge structures. The main challenges arising in this context relate to the high dimensionality of the damage condition space, which will require more efficient distributional models such as copulas.

Appendix A: Exploring results when loading conditions are unavailable

This appendix considers the case when the loading condition features describing the wind and wave excitation are unavailable. We train our inverse model ℐ_θ, ignoring these features for the case where only the roll DOF is available. The design and training specifications are similar to those summarized in Table 1. Figure A1 presents the results obtained for the first four test case scenarios. Analyzing the results, we observe that, although the distributional model partially captures the damage condition in some cases (see, e.g., Fig. A1a, c, and g), in other scenarios it is far from producing reliable outcomes. The lack of information regarding the loading conditions misleads the condition estimates, mainly for variable z₁ (anchoring). For example, in Fig. A1b and e, we observe that the distributional model spreads along the anchoring axis (z₁) with no identification of the existing multimodality.

Figure A1 First four test examples for 1 DOF when omitting wind and wave condition measurements in the inverse training stage.

Appendix B: Exploring results for simultaneous damage in two mooring lines

In this appendix, we analyze the detection capability of simultaneous damage occurring in two different mooring lines. To remain in the two-dimensional space, we consider only anchoring damage. In particular, this section explores anchoring damage to a lateral mooring line, as well as the downstream line parallel to the surge motion. The parameterization of anchoring degradation for the platform's lateral line remains unchanged with respect to the main body of this study. Anchoring damage displaces the anchor of the floater's downstream line up to 20 m further in the direction of the wind–wave excitation. Figure B1 presents the probability and data misfit distributions for four test cases, considering as input data the 3 rotation DOFs (roll, pitch, and yaw). The figure shows the robust performance of the model to concurrent damage in various mooring lines. However, if we constrain the input data to the roll DOF only, we find that the measured features are unable to identify the damage in the second mooring line. Figure B2 illustrates this for one test scenario. The results highlight the importance of selecting damage-sensitive features to capture the different modes of failure that may occur in the system.

Figure B1 Four test cases for simultaneous anchoring damage on two mooring lines considering roll, pitch, and yaw motion. The left-hand figures represent the contour plot of the estimated probability density functions $q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$. The right-hand figures represent the data misfit (ℒ_ℳ) value.

Figure B2 Test case for simultaneous anchoring damage on two mooring lines considering only roll motion.

Appendix C: Exploring the generalization capability for unseen-damage scenarios

This appendix explores the potential of our proposed method in evaluating damage scenarios different from the training cases. Given the extrapolation limitations of neural networks, we study scenarios of the same nature (the same type of damage occurs, namely anchoring and biofouling), although with a remarkable reduction in severity. We perform two different analyses for the inverse: (i) training only with scenarios where damage z₁ (anchoring) is always ≥20 % and testing with scenarios with z₁<20 % and (ii) training only with scenarios where damage z₂ (biofouling) is always ≥20 % and testing with scenarios with z₂<20 %. In both cases, we consider using the 3 rotation DOFs as the input measured signals (roll, pitch, and yaw).

Figure C1 shows four randomly selected test cases for the first analysis (anchoring kept ≥20 % during training). We observe that the distributional model successfully approximates the true solution, with the uncertainty spreading mostly in the direction of the constrained feature (i.e., anchoring) within acceptable thresholds. Hence, the developed method can inform the damage condition and indicate the direction of uncertainty, even for cases where the evaluated scenario lies outside the training domain.

Analogously, Fig. C2 shows four test scenarios for the second analysis (biofouling kept ≥20 % during training). The results reveal that the highest density values concentrate close to the true solution, although with some error (see, e.g., Fig. C2c, where the highest PDF values are shifted towards smaller anchoring damage, but the level of biofouling damage is well captured).

In summary, the results demonstrate that despite these not being included during the training phase, the DNN successfully identifies damage scenarios with lighter severity. The analyses demonstrate the generalizability of the method for damage scenarios that have the same nature but are significantly different in severity from those considered during training, as well as the potential of the distributional model to indicate the direction of the uncertainty.

Figure C1 Posterior PDF ($q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$) contour plots for four test scenarios when training with anchoring (z₁) scenarios ≥20 %.

Figure C2 Posterior PDF ($q_{\mathit{\theta }}\left(\mathbf{z}\right|\mathbf{m},\mathbf{w})$) contour plots for four test scenarios when training with biofouling (z₂) scenarios ≥20 %.