Measuring the Robustness of Optimal Design Solutions for Wave Energy Converters via a Stochastic Approach

Abstract

Among the challenges generated by the global climate crisis, a significant concern is the constant increase in energy demand. This leads to the need to ensure that any novel energy systems are not only renewable but also reliable in their performance. A viable solution to increase the available renewable energy mix involves tapping into the potential available in ocean waves and harvesting it via so-called wave energy converters (WECs). In this context, a relevant engineering problem relates to finding WEC design solutions that are not only optimal in terms of energy extraction but also exhibit robust behavior in spite of the harsh marine environment. Indeed, the vast majority of design optimization studies available in the state-of-the-art consider only perfect knowledge of nominal (idealized) conditions, neglecting the impact of uncertainties. This study aims to investigate the information that different robustness metrics can provide to designers regarding optimal WEC design solutions under uncertainty. The applied methodology is based on stochastic uncertainty propagation via a Monte Carlo simulation, exploiting a meta-model to reduce the computational burden. The analysis is conducted over a dataset obtained with a genetic algorithm-based optimization process for nominal WEC design. The results reveal a significant deviation in terms of robustness between the nominal Pareto set and those generated by setting different thresholds for robustness metrics, as well as between devices belonging to the same nominal Pareto frontier. This study elucidates the intrinsic need for incorporating robust optimization processes in WEC design.

1. Introduction

In 2020, the United Nations emphasized the urgent need for every nation to address a critical goal for the near future: achieving carbon neutrality by 2050 [1]. However, the pursuit of this objective faces a significant obstacle due to the rising energy demand, which is currently largely reliant on fossil fuels [2]. For many years, various renewable energy sources have been under the scope of bridging this gap, with ocean wave energy emerging as a large and almost untapped resource. The global wave energy resource is estimated to lie in the range of 1 to 10 TW, but the exact extractable wave power remains the subject of ongoing research. Theoretical calculations, as presented in [3], suggest that this global resource totals around 32,000 TWh/y with a mean power output of $3.65$ TW. Using the same assessment methodology and data, the Oceanic Energy System and the International Renewable Energy Agency report a value of 29,500 TWh/y [4,5]. These numbers demonstrate that wave energy could play a pivotal role in the race towards achieving carbon neutrality by 2050.

From a historical perspective, the origins of technology designed to harness wave energy date back to 1799 in France.1 Today, wave energy is gaining attention due to its potential to become a competitive new player in the renewable energy market. One of the primary reasons for this is the absence of marketable devices for energy production, essential for establishing viable wave energy, mainly due to the lack of conversion concepts’ convergence. A broad overview of the topic can be found in [6].

Wave energy presents some advantages (therefore sparking interest) when its power density is compared to that of other well-known renewable energy sources. Waves are mainly generated via the interaction between wind and the sea surface, in which the constant mechanical action of the wind, acting as a tangential effort, leads to the formation of waves. On the other hand, the deployment of wave energy converter (WEC) technology has yet not reached a competitive level compared to other clean energy industries. Different key points have been mapped to proceed towards WEC’s successful development and commercialization. Concerning the latter topic, various studies (see, e.g., [6,7]) underline that, for a WEC investment to be considered appealing to stakeholders and investors, it should meet the following criteria: it should offer energy production costs that are competitive on the market; it should present a secure and dependable investment opportunity; it should demonstrate reliability in supporting grid operations; it should contribute positively to society or have societal benefits; it should align with permitting and certification requirements; it should prioritize operational safety; and it should have the potential for global deployment. A general aspect that encompasses all the issues listed immediately above is the common emphasis on the techno-economic reliability of WEC technology. This, in fact, is consistent with the conclusions argued for in [6]: to reduce investor risk, it is necessary to reduce both the levelized cost of energy (LCoE) and its associated uncertainty. The highlighted challenges are classically addressed via robust optimization (RO). RO frameworks aim to find solutions to optimization problems while incorporating uncertainty in the evaluation of the fitness function. RO techniques’ purpose is to achieve optimal solutions that are less prone to being affected by uncertainty.

To illustrate this, in Figure 1, it is possible to distinguish among different local optima, i.e., the minima of the objective function $f{\left(x\right)}_{obj}$. Each one is characterized by a different behavior. Point B is referred to as a robust optimum since it remains locally optimal even when a relatively large variation $\Delta x$ occurs. Conversely, despite the fact that A is a global optimum, if the same variation $\Delta x$ perturbs point A, the objective function increases sharply, losing optimality with respect to point B; i.e., it is effectively more sensitive to variations in x. A conventional formulation of the associated RO problem is provided in Section 2.

The developments presented within this study were effectively motivated by the limited literature available on the analysis of the performance robustness of WECs from a design perspective. In particular, and with respect to the classification of robustness metrics outlined in a previous paper [9], it can be asserted that there exists a lack of comprehensive and systematic research in the field of WEC design addressing all four subcategories that underlie the idea of robustness. The first subcategory is the so-called robustness of a concept, and it relies on robustness metrics that assess the impact of design parameters or independent noise factors on the investigated performance. Subsequently, when focusing on the size of the feasible design space, the robustness of a design and its indicators are defined. Another branch of robustness measures assesses the different expectations and dispersion of the system’s output, i.e., the robustness of a function. Finally, the robustness of a product classification exploits metrics that assess the likelihood of satisfactory performance conformity under the influence of ongoing variations. To date, within the realm of wave energy conversion, the literature on robust optimization has predominantly focused on the development of robust control frameworks [10,11,12,13,14,15,16,17]. Furthermore, concerning the study of wave energy and WEC uncertainties, the existing literature primarily covers environmental [18,19,20,21,22] and full/high-scale prototyping or experimental uncertainties and subsystems’ influence on performance [23,24,25,26,27,28,29]. The neighboring Reliability-Based Design Optimization (RBDO) field of research directs its attention to WEC structural and maintenance cost uncertainties [30] (in which the so-called Wavestar device is considered as a case study) and the power take-off (PTO) reliability relationship with hull geometry [31]. A review of RBDO studies and their application in the offshore renewable energy sector is presented in [32].

Motivated by the necessity of WEC systems’ RO, the present study evaluates different robustness indices (RIs) for WEC design solutions and assesses their impact on devices’ overall normalized power performance (Equation (18); i.e., it performs a stochastic Robustness Quantification (RQ) investigation for WEC systems. To achieve this objective, the present study conducted an analysis to gather information regarding the robustness of these devices generated in terms of a dataset processed during a multi-objective nominal optimization process. First, nominal Pareto set devices RI were evaluated, and information about their statistical behavior was obtained accordingly. This nominal classification is discussed in accordance with the standard Pareto criterion. Threshold criteria were then imposed on robustness indices and new Pareto research, highlighting the main differences between obtained and nominal limit sets. Furthermore, due to the computational burden associated with the uncertainty propagation, performed via a Monte Carlo simulation, a surrogate-assisted RQ method is proposed. The reason for presenting the results and devices’ data for the case study examined in normalized form is the confidential and industrial nature of the analyzed technology. Despite this limitation, the results presented in this non-dimensional form are sufficient to draw conclusions regarding robustness and, thus, adequate for the purpose of this study.

The remainder of the paper is structured as follows. In Section 2, an overview of the WEC robust optimization problem is provided, together with a description of different RQ approaches and robustness indices. Section 3.2 and Section 3.1, respectively, discuss and provide details on the specific framework applied in this study and the applicative case of analysis, i.e., the Inertial Sea Wave Energy Converter (ISWEC) device. Section 4 presents the main results and discussion. Finally, in Section 5, key conclusions of this study are provided, including future lines of research.

2. Robustness Quantification

In this section, a brief overview of robust optimization and its formulation is presented, together with a brief review of different robust quantification metrics.

2.1. Robust Optimization Overview

In the literature, the robustness concept has been articulated in different ways, according to different fields of study. One of the first definitions of robustness was provided by Taguchi for the quality improvement of industrial processes: “a product whose performance is minimally sensitive to factors causing variability (at the lowest possible cost)” [33]. The robustness concept also plays a central role in the already-mentioned robust control field, where different applications can be found. In this specific case, the aim is to design and synthesize a control a strategy that is ideally insensitive to either parametric or non-parametric uncertainty, able to guarantee both a prescribed set of performance specifications (e.g., maximizing energy harvesting) and the robust stability of the controlled system [34]. Regarding WEC control, in [10], a hierarchical robust strategy was devised to mitigate the controller susceptibility to modeling inaccuracies and potentially unmodeled nonlinear influences. Model sensitivity issues for different WEC control system architectures are addressed in [11], leveraging the inherent parametric uncertainty in standard device modeling. Moreover, a moment-based framework to design robust energy-maximizing optimal controllers for WEC is presented in [12].

The pivotal role of robustness is also clear in the context of optimization problems, in which the RO research field is relevant. As a historical note, consider that in 1973, Soyster [35] formulated a linear programming (LP) optimization model, yielding feasible solutions for all data belonging to a convex set. However, this approach leaned towards excessive conservatism, compromising optimality for enhanced robustness. In subsequent studies by Ben-Tal, Nemirovski, El-Ghaoui, and Lebert [34,36], spanning from 1997 to 2000, less conservative solutions were proposed. These strategies tackled uncertain linear problems with ellipsoidal uncertainties, requiring the solution of robust counterparts in a conic quadratic form. Despite offering reduced conservatism, this approach introduced practical drawbacks, as it led to inherently nonlinear (yet convex) formulations with increased computational demand compared to Soyster’s earlier linear framework. Expanding on this historical note, Bertsimas and Sim [37] integrated the insights from Soyster and Ben-Tal et al. into a fully controllable framework. This latter approach allowed for the adjustment of conservatism levels with probabilistic guarantees. Notably, their proposed robust counterparts were linear optimization problems, enabling a seamless generalization to discrete optimization problems. A key novelty was the introduction of the so-called $\mathsf{\Gamma }$ parameter, which refers to a protective threshold that plays a relevant role in fine-tuning robustness against conservatism levels.

The second main methodology for solving RO problems is called “Robust Design Optimization” (RDO) [38]. This method, complementary to RO, employs a stochastic approach based on different robustness metrics, aiming to achieve a system design insensitive to uncertainty, external noises, perturbations, potential model inaccuracies, and design tolerances [8]. The latter is of significant importance within the realm of simulation-based optimization, particularly when addressing complex engineering challenges, such as the design of WECs. Typically, the RDO framework is based on the RQ process, which is a stochastic approach to uncertainty propagation that can be outlined in the following main steps.

• Parameter uncertainty vectors’ ($\overrightarrow{u}$) and design ($\overrightarrow{x}$) parameter vectors’ definition:
  • $\overrightarrow{u}={[u_1\dots u_N]}^{⊺}\in {\mathbb{R}}^N$;
  • $\overrightarrow{x}={[x_1\dots x_D]}^{⊺}\in {\mathbb{R}}^D$,
  where N is the number of system parameters affected by uncertainty, whose magnitude is described by the relative $n^{th}$-element of the uncertainty vector $\overrightarrow{u}$, and D is the dimension of the design parameter vector $\overrightarrow{x}$, i.e., the number of decision elements that define each possible optimization problem solution.
• Uncertainty propagation via sampling techniques.
• Robustness metric quantification.

A description of the proposed RQ framework for the present study is provided in Section 3.2. Therefore, it is common that, in the RDO field, the objective function is not known exactly but is evaluated in a computationally expensive fashion, leveraging stochastic logic. Consequently, efforts have been directed towards studies concerning evolutionary algorithms [39,40], especially via multi-objective optimization, to find a robust Pareto frontier [41,42].

2.2. Problem Formulation and Robustness Indices

As discussed within Section 1, this paper is focused on the realm of RDO, and it is primarily aimed at serving as a quantitative exploration of the challenges involved when incorporating uncertainty in optimal WEC design, together with any consequent implications. In particular, a stochastic approach was adopted, and the motivation behind this choice lay in avoiding the adoption of assumptions potentially restrictive to the system and exploring computer simulation-experiments’ viability. Therefore, a general RDO problem formulation can be mathematically stated as follows:

$$\underset{\overrightarrow{x}\in \mathcal{X}}{min}\overrightarrow{f}(\overrightarrow{x},\overrightarrow{u}),$$

(1)

which is subject to the following:

$$\begin{array}{c}h(\overrightarrow{x},\overrightarrow{u})=0,\\ g(\overrightarrow{x},\overrightarrow{u})<0,\end{array}$$

(2)

where the following points apply:

• ${\overrightarrow{u}}^{min}<\overrightarrow{u}<{\overrightarrow{u}}^{max}$ are the upper and lower bounds of the uncertainty parameters’ vector. The space resulting from limiting $\overrightarrow{u}$ can be referred to as the parameter uncertainty space $\mathcal{U}$.
• ${\overrightarrow{x}}^{min}<\overrightarrow{x}<{\overrightarrow{x}}^{max}$ are the upper and lower bounds of the design vector. The space resulting from limiting $\overrightarrow{x}$ can be referred to as the parameter design space $\mathcal{X}$.
• $h(\overrightarrow{x},\overrightarrow{u})\in {\mathbb{R}}^q$ represent q equality constraints.
• $g(\overrightarrow{x},\overrightarrow{u})\in {\mathbb{R}}^k$ represent k inequality constraints.
• $\overrightarrow{f}(\overrightarrow{x},\overrightarrow{u})={[f_1(\overrightarrow{x},\overrightarrow{u})\dots f_m(\overrightarrow{x},\overrightarrow{u})]}^{⊺}\in {\mathbb{R}}^m$ is the vector containing the potentially different m objective functions.

A comprehensive survey, centered on the methods for solving the RDO problem and the associated robustness quantification metrics, is presented in [38]. In particular, RDO problems are solved via single-objective robust design optimization or multi-objective robust design optimization. Single-objective RDO frameworks involve setting the robustness metric as part of the optimization objective function. On the other hand, multi-objective RDO processes incorporate the RI as an additional objective function [43]. Furthermore, another option is to rely on supplementary constraints to enhance robustness, i.e., RBDO [44]. A well-chosen robustness measure offers a valuable means to discern a system’s responsiveness (or lack thereof) to factors that might result in sub-optimal performance. A glimpse of some of the robustness metrics employed in the literature [43,45], and considered within the present study, is provided in

Table 1. Robustness metrics and indices.

Definition	Index
(I) Average value of the $m^{th}$-performance probability density function (PDF).	${\widehat{\mu }}_i=\frac{1}{N_{sample}}{\sum }_{p=1}^{N_{sample}}\frac{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}=\frac{{\mu }_i(\overrightarrow{x_i},\overrightarrow{u})}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}$    , where $N_{sample}$ is the quantity of samples used in order to evaluate the discrete PDF mean value. The same formulation results are valid for all the other Definitions in the present table.
(II) Standard deviation of the $m^{th}$-performance PDF.	${\widehat{\sigma }}_i=\frac{1}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}\sqrt{\frac{{\sum }_{p=1}^{N_{sample}}{[f_m(\overrightarrow{x_i},\overrightarrow{u_p})-{\mu }_i(\overrightarrow{x_i},\overrightarrow{u})]}^2}{N_{sample}-1}}$.
(III) The “k-sigma” approach.	${\widehat{\mu }}_i+k{\widehat{\sigma }}_i$, where k is an integer number commonly set equal to 3 or 62
(IV) Symmetric robustness index.	$R_i=\frac{s_{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}+\left|{\mu }_{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}-f_m(\overrightarrow{x_i},\overrightarrow{u_0})\right|}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}$, with: $s_{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}=\sqrt{\frac{{\sum }_{p=1}^{N_{sample}}{[f(\overrightarrow{x_i},\overrightarrow{u_p})-f_m(\overrightarrow{x_i},\overrightarrow{u_0})]}^2}{N_{sample}-1}}$ and ${\mu }_{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}=\frac{1}{N_{sample}}{\sum }_{p=1}^{N_{sample}}{f}_m(\overrightarrow{x_i},\overrightarrow{u_p})$.
(V) Asymmetric robustness index.	$Q_i=\frac{Q_0}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}$.
(VI) The upper limit of the $m^{th}$-performance PDF.	$sup{\widehat{f}}_i=\underset{\overrightarrow{u}}{max}\frac{f_m(\overrightarrow{x_i},\overrightarrow{u_p})}{f_m(\overrightarrow{x_i},\overrightarrow{u_0})}\forall p\in [1\dots N_{sample}]$.

, together with their formulations. The notation used in Table 1 refers to the so-called zero-uncertainty vector $\overrightarrow{u_0}$ as the parameters uncertainty vector with all elements equal to zero, i.e., $\overrightarrow{u_0}={[0\dots 0]}^{⊺}\in {\mathbb{R}}^N$. Therefore, the performance evaluated for $\overrightarrow{x_i}$ and $\overrightarrow{u_0}$ ($f_m(\overrightarrow{x_i},\overrightarrow{u_0})$) is named the $m^{th}$-performance nominal expected value.

Definition I in Table 1 is consistent with the idea of optimizing the average $m^{th}$-performance of the analyzed system. This option could direct the final choice (as well as the process of exploration of the optimization itself) towards individuals characterized by good average performance but potentially highly dispersed around the nominal expected value.

Following Definition II in Table 1 means considering only the samples’ distribution concerning the mean value of the PDF. Due to its limited possibility to describe a system’s performance in terms of an absolute value, ${\widehat{\sigma }}_i$ is mainly adopted in multi-objective RDO frameworks as part of the objective functions, together with other parameters, e.g., ${\widehat{\mu }}_i$, or the performance nominal value $f(\overrightarrow{x_i},\overrightarrow{u_p})$ directly.

The “k-sigma” approach establishes a range of performance deviations from the PDF’s statistical average value by merging the approaches specified in Definitions I and II within Table 1. This is accomplished by applying a multiplicative coefficient on ${\widehat{\sigma }}_i$. The smaller this deviation, the more robust the individual examined.

In Definition IV (Table 1), $R_i$ is defined as “symmetric” since its formulation weights every possible deviation from the nominal expected value. In particular, a smaller value of $R_i$ indicates a higher level of robustness in the system being considered. The framework in which $R_i$ is usually used can involve its optimization, together with the robust performance objective. Thanks to the influence of both the distance between the average value of the distribution and its spread with respect to the expected nominal value $f_m(\overrightarrow{x_i},\overrightarrow{u_0})$, Definition IV is well suited for metrics for which the purpose of the designer is to keep the performance very close to a specific target.

On the other hand, the $Q_i$ “asymmetric” robustness index (Definition V in Table 1) penalizes only the deviations from the nominal value that result in poorer performances and, consequently, worse values of the objective function in question. In the case in which the $m^{th}$ performance needs to be minimized, $Q_0$ represents the objective function value below the $q\%$ (typically greater than $90\%$) of observed occurrences. In this scenario, a smaller value of $Q_i$ indicates a better risk measure. Conversely, if the $m^{th}$ performance needs to be maximized, $Q_0$ is the objective function value above the $q\%$ of observed occurrences. In this case, a larger value of $Q_i$ indicates a better risk measure. To enable a comparison across different devices, the nominal value $f_m(\overrightarrow{x_i},\overrightarrow{u_0})$ is used in the index for normalization purposes. The $Q_i$ index offers insights into device quality concerning the $m^{th}$ objective function. Given this capability, it can be used for systems belonging to various stages of product development and their potential for exploitation. Nonetheless, this measurement does not inform the designer about specific occurrences’ values; it reports only cumulative information.

To perform a robustness quantification that considers the upper limit in a minimization problem (respectively, the lower limit in a maximization case), Definition VI in Table 1 can be adopted, i.e., the Worst-Case Scenario (WCS). The WCS RI presents two main drawbacks: over-conservative solutions are often reached, and an extensive exploration phase of the uncertain parameter’s N-dimensional space (to evaluate the real WCS) is required, leading to computationally expensive evaluations. The interested reader can find comprehensive examinations of robustness metrics employed in RDO processes in the works of Beyer et al. [38] and Moritz Göhler et al. [9]. Moreover, in [41,42], Deb et al. define two different robust solution approaches to optimization, focusing on evolutionary algorithms.

When quantifying robustness, it is crucial to take into account the nature of uncertainties. Typically, in the realm of RDO, the primary types of uncertainties under scrutiny are epistemic and aleatory. The first type of uncertainty, i.e., epistemic uncertainty, stems from a lack of knowledge that can be mitigated by acquiring more information about the parameter. For instance, in the context of renewable energy systems, epistemic uncertainty in the typical annual electricity demand may arise from unknown occupant behavior. On the other hand, aleatory uncertainty is associated with the unpredictable variation of parameter values and is, therefore, inherent and irreducible. For example, aleatory uncertainty regarding wholesale electricity prices results from the unpredictable variations of this parameter in the future. Both types of uncertainties can be characterized using interval-valued statistical moments [47].

It is readily evident that adopting a stochastic approach to propagating uncertainties in quantifying the robustness of the previously described system presents several challenges, with the most prominent being elevated computational expenses. To address this challenge, so-called surrogate-assisted RDO has been developed. These methods rely on the use of mathematical models that approximate a time-demanding model to reduce the computational cost of the objective function (i.e., a so-called surrogate model or metamodel), making the computational expense required to solve the optimization problem tractably. Various types of mathematical models are employed in the literature, with some of the most common ones including Gaussian process regression (GPR), the Kriging Model, and polynomial chaos expansion [48,49]. An additional avenue for reducing the computational cost of RDO processes could involve a shift in the problem-solving approach, transitioning from a stochastic framework to a deterministic framework. Assuming the system exhibits behavior such that the WCS relies on the extremes of the uncertainty-affected parameter range, the designer could choose to evaluate the objective function solely at the vertices of the resulting N-dimensional polytope. For a single uncertain variable, there would be two vertex points (the uncertain space will be described by a single line); for two variables, there would be four evaluation points (forming a square); for three variables, eight points (describing a cube), and so forth. This approach would result in a reduction in computational costs by diminishing the number of expensive evaluations required. The advantages and disadvantages associated with the various approaches outlined earlier are condensed in

Table 2. Main advantages and disadvantages of different approaches to optimization in the presence of uncertainty, inspired by [50].

Method	Advantage	Disadvantage
Stochastic RDO	Provide trade-off solutions and detailed information on the output under examination distribution.	Very computationally expensive.
Deterministic RDO	Less computationally expensive.	Need the assumption of good behavior in the N-dimensional polytope. Overly conservative.

.

3. Methodology

The ISWEC device and its model under analysis are described in more detail in the present section, along with a technical introduction to the technology. Furthermore, the objectives of the analyses conducted, the framework, and the methodologies used to achieve them are introduced. The main goal of this work was to compare different stochastic robustness metrics and evaluate the information that each of them provides to the designer. The robustness metrics used included the six previously introduced in Table 1, namely ${\widehat{\mu }}_i$, ${\widehat{\sigma }}_i$, $sup\left({\widehat{f}}_i\right)$, ${\widehat{\mu }}_i+k{\widehat{\sigma }}_i$, $R_i$ and $Q_i$.

3.1. Model

The ISWEC device (Figure 2) is a WEC characterized by a pitch-floating hull, featuring a fully enclosed design. The device was designed to perform as well as the Pendulum Wave Energy Converter device [51] in the demanding conditions of the Mediterranean Sea, and its performance was optimized with this aim. A concise and general overview of ISWEC and its operational principles is provided below.

The heart of this system lies within the floater, housing a PTO mechanism that leverages the precession motion of a gyroscope system to capture the mechanical energy induced via wave action. This precession motion is initiated through the dynamic coupling between pitch motion and flywheel rotation. To prevent gyroscope stabilization, an eccentric mass is integrated, serving the purpose of inducing elastic recoil. For the enhancement of energy harvesting and resonance tuning, a well-tailored PTO control logic can be incorporated [52]. In the scope of this study, the PTO control logic description is provided in subsequent paragraphs. In the pursuit of investigating the role of uncertainty in the WEC design process, it is imperative to conduct a series of numerical experiments through simulations. Consequently, consideration must be given to the selection of an appropriate numerical model. A state-of-the-art perspective on the topic of numerical modeling applied to the WEC case is presented in [53].

In this section, the ISWEC modeling details are described. In Figure 3, the ISWEC functional diagram is presented. The block chart provides a concise visual representation of the device’s behavior, encompassing forces and reactions between the floater and the gyroscope, as well as the control actions and signal flow of the PTO system [54]. In the image, we can recognize the following elements:

• ${\overrightarrow{f}}_f$ is the combined gravitational and fluid forces acting on the floater.
• $T_{\epsilon }$ and ${\overrightarrow{f}}_{gf}$ represent the gyroscopic effects exchanged between the floater and the gyroscope units.
• $T_{PTO}$ is used to represent the PTO reaction acting on the gyroscope shaft.
• ${\mathbf{X}}_f$, ${\mathbf{X}}_{\epsilon }$, and ${\mathbf{X}}_{PTO}$ describe the motion of the floater, gyroscope, and PTO, respectively.
• ${\mathbf{C}}_{PTO}$ stands for the control signals generated via the ISWEC’s control system.

The mathematical equations that govern the ISWEC model were derived from the coupling of the WEC hydrodynamic with the gyroscope and PTO dynamics. As suggested in [54], a classic time domain representation is well suited for a more realistic and advanced simulation phase of WEC development than the frequency domain approach; that is, via the time domain, it is possible to consider non-linearities in the model [55]. Following the formulation presented in [54], the time domain system was composed with the involvement of the floater equation and the radiation forces’ approximation, together with the gyro-pendulum model and the consideration of gyroscopic effects. Therefore, the non-linear state-space Time Domain Model (TDM) equations can be written in the following form:

$${\mathbf{M}}_T{\ddot{\mathbf{X}}}_T\left(t\right)+{\mathbf{B}}_T{\dot{\mathbf{X}}}_T\left(t\right)+{\mathbf{K}}_T{\mathbf{X}}_T\left(t\right)+{\mathsf{\Theta }}_T({\ddot{\mathbf{X}}}_T\left(t\right),{\dot{\mathbf{X}}}_T\left(t\right),{\mathbf{X}}_T\left(t\right))={\mathbf{f}}_T\left(t\right).$$

(3)

For a system of dimension $n_T$, the TDM model formulation of Equation (3) includes the kinematic variable vector ${\mathbf{X}}_T\left(t\right)\in {\mathbb{R}}^{n_T}$, the mass matrix ${\mathbf{M}}_T\in {\mathbb{R}}^{n_T\times n_T}$, the damping matrix ${\mathbf{B}}_T\in {\mathbb{R}}^{n_T\times n_T}$, the stifness matrix ${\mathbf{K}}_T\in {\mathbb{R}}^{n_T\times n_T}$, the non-linear function ${\mathsf{\Theta }}_T:{\mathbb{R}}^{n_T}\to {\mathbb{R}}^{n_T}$, and the external time-dependent excitation forces ${\mathbf{f}}_T\left(t\right)\in {\mathbb{R}}^{n_T}$. Several non-linearities affect the term ${\mathsf{\Theta }}_T$, deriving from the floater, gyro-pendulum and PTO subsystems.

Alternatively, Frequency Domain Models (FDMs) offer a computationally efficient solution via the boundary element method. To adopt an FDM means to accept some assumptions related to the realistic non-linear system behavior. The two presupposition characteristics of the ISWEC FDM regard the gyroscopic torques and pendulum elastic recall and the exclusion of saturations from the PTO model. That is, the ISWEC system is assumed to work only in the linear PTO model zone and the gyroscopic torques and pendulum elastic recall are linearized around the equilibrium point of the precession axis. After those assumptions, the model can be written in a matrix form concerning a harmonic input at the intended angular frequency $\omega$:

$$[-{\omega }^2{\mathbf{M}}_F\left(\omega \right)+i\omega {\mathbf{B}}_F\left(\omega \right)+{\mathbf{K}}_F]{\widehat{\mathbf{X}}}_F={\mathbf{H}}_{\eta F}\left(\omega \right)\widehat{\eta },$$

(4)

where ${\mathbf{M}}_F\left(\omega \right)\in {\mathbb{R}}^{n_F\times n_F}$ is the frequency-dependent mass matrix (which includes the floater mass matrix, the gyroscope inertia, and the added mass contribution), ${\mathbf{B}}_F\left(\omega \right)\in {\mathbb{R}}^{n_F\times n_F}$ is the frequency-dependent damping matrix, ${\mathbf{K}}_F\in {\mathbb{R}}^{n_F\times n_F}$ is the stiffness matrix, and ${\mathbf{H}}_{\eta F}\left(\omega \right)\in {\mathbb{C}}^{n_F}$ is the complex transfer function between the wave profile $\eta$ and the resultant excitation forces F, while $n_F$ is the dimension of the frequency domain system. The FDM presentation also allows us to fully describe the statistical properties of ISWEC outputs in terms of the power spectral density (PSD) of WEC displacements ${\mathbf{S}}_{X_F{X}_F}\left(\omega \right)$ [54]. The latter statement is explained because, when a Gaussian stochastic input passes through a linear transfer function, the output also follows a Gaussian distribution [53]. Furthermore, the PDF of sea wave elevation $\eta$ can be considered to be in accordance with the assumption of a Gaussian-distribution sea [56]. Therefore, when adopting the FDM representation, Gaussian PDF properties can be assumed for both floater and gyroscope displacements and evaluated, via the transfer function ${\mathbf{H}}_{\eta {\mathbf{X}}_F}\left(\omega \right)$, using Equation (6) [54]:

$${\mathbf{S}}_{X_F{X}_F}\left(\omega \right)={\mathbf{H}}_{\eta X_F}\left(\omega \right)S_{\eta \eta }\left(\omega \right){\mathbf{H}}_{\eta X_F}^{\ast }\left(\omega \right).$$

(5)

In Equation (5), the notation $\ast$ stands for the complex conjugate operator. Moreover, $S_{\eta \eta }$ represents the wave elevation PSD.

$${\mathbf{H}}_{\eta X_F}\left(\omega \right)=\frac{{\mathbf{H}}_{\eta F}\left(\omega \right)}{-{\omega }^2{\mathbf{M}}_F\left(\omega \right)+i\omega {\mathbf{B}}_F\left(\omega \right)+{\mathbf{K}}_F}.$$

(6)

In order to incorporate non-linear factors to enhance the model’s accuracy, the Spectral Domain Model (SDM) representation can be considered. SDM’s purpose is to obtain an approximate solution of a generic TDM while avoiding the extensive computational overhead associated with it. SDM’s practical implementation is accomplished by describing the system as probabilistic, considering its inputs as a stochastic ergodic process [53]. Hence, as observed earlier in the context of an FDM, when the wave spectrum is subjected to a suitable transfer function, it generates a probabilistic prediction of the WEC response [54]. Therefore, in terms of the PSD of WEC displacements, the following applies:

$${\mathbf{S}}_{X_S{X}_S}\left(\omega \right)={\mathbf{H}}_{\eta X_S}\left(\omega \right)S_{\eta \eta }\left(\omega \right){\mathbf{H}}_{\eta X_S}^{\ast }\left(\omega \right),$$

(7)

where ${\mathbf{H}}_{\eta X_S}\left(\omega \right)\in {\mathbb{C}}^{n_S}$ is the SDM transfer function between WEC displacement and wave elevation, $n_S$ is the dimension of the spectral domain system, and ${\mathbf{S}}_{X_S{X}_S}\left(\omega \right)$ is the ISWEC outputs’ PSD. In SDM, the generic system equations can be written as follows:

$${\mathbf{M}}_S{\ddot{\mathbf{X}}}_T\left(t\right)+{\mathbf{B}}_S{\dot{\mathbf{X}}}_T\left(t\right)+{\mathbf{K}}_S{\mathbf{X}}_T\left(t\right)={\mathbf{f}}_T\left(t\right),$$

(8)

where the ${\mathbf{X}}_T$ vector identifies the hull’s motion (${\mathbf{X}}_f$), and the matrix ${\mathbf{M}}_S\in {\mathbb{R}}^{n_S\times n_S}$ is the equivalent mass matrix, ${\mathbf{B}}_S\in {\mathbb{R}}^{n_S\times n_S}$, is the equivalent damping matrix, and ${\mathbf{K}}_S\in {\mathbb{R}}^{n_S\times n_S}$ is the equivalent stiffness matrix.

Moreover, the following applies:

$$\begin{array}{c}{\mathbf{M}}_S={\mathbf{M}}_T+{\mathbf{M}}_{eq},\end{array}$$

(9)

$$\begin{array}{c}{\mathbf{B}}_S={\mathbf{B}}_T+{\mathbf{B}}_{eq},\end{array}$$

(10)

$$\begin{array}{c}{\mathbf{K}}_S={\mathbf{K}}_T+{\mathbf{K}}_{eq}.\end{array}$$

(11)

These matrices are composed of a first linear part related to the linear terms of Equation (3) and a second part that represents the statistical time-varying behavior of the time domain system’s non-linear terms. Therefore, this second part of each matrix is defined using the gradient and the expected value $\langle .\rangle$ operator as follows:

$$\begin{array}{c}{\mathbf{M}}_{eq}=\langle {\nabla }_{\ddot{\mathbf{X}}}\mathsf{\Theta }\rangle ,\end{array}$$

(12)

$$\begin{array}{c}{\mathbf{B}}_{eq}=\langle {\nabla }_{\dot{\mathbf{X}}}\mathsf{\Theta }\rangle ,\end{array}$$

(13)

$$\begin{array}{c}{\mathbf{K}}_{eq}=\langle {\nabla }_{\mathbf{X}}\mathsf{\Theta }\rangle .\end{array}$$

(14)

where the ${\nabla }_{\mathbf{X}}$ operator defines the gradient with respect to $\mathbf{X}$. The SDM approach involves the minimization of the mathematical expectation of the discrepancy between the non-linear system and its linearized counterpart, which is represented as $e_S$:

$$e_S={\mathsf{\Theta }}_T({\mathbf{X}}_T\left(t\right),{\dot{\mathbf{X}}}_T\left(t\right),{\ddot{\mathbf{X}}}_T\left(t\right))-{\mathbf{M}}_{eq}{\ddot{\mathbf{X}}}_T\left(t\right)-{\mathbf{B}}_{eq}{\dot{\mathbf{X}}}_T\left(t\right)-{\mathbf{K}}_{eq}{\mathbf{X}}_T\left(t\right).$$

(15)

After the statistical linearization of the system is performed, it is possible to evaluate the SDM transfer function:

$${\mathbf{H}}_{\eta X_S}\left(\omega \right)=\frac{{\mathbf{H}}_{\eta S}\left(\omega \right)}{-{\omega }^2{\mathbf{M}}_S\left(\omega \right)+i\omega {\mathbf{B}}_S\left(\omega \right)+{\mathbf{K}}_S}.$$

(16)

This problem lacks an analytical solution, necessitating an iterative procedure, which is presented in [54]. The use of the spectral model is justified as long as the assumption of the system’s Gaussian behavior remains valid. If the previously described assumption were to fail (e.g., in the presence of dominant non-linear forces), the system would exhibit non-Gaussian behavior, making the spectral model inaccurate and erroneous.

While both spectral and frequency domain models may sacrifice some degree of accuracy compared to time domain simulations, they are recommended for analyses seeking to compare different WEC devices and design solutions, as suggested in previous works [54] and the literature [53]. For the forthcoming analyses, we opted for the spectral domain as our mathematical model of choice. Given the importance of the process of converting wave energy into the electrical energy of the PTO, a brief description of its model is provided below. The mechanical PTO consists of a one-stage gearbox connecting the gyroscope shaft to the electrical generator. The control logic used in this work for the PTO was founded on the principle of impedance matching, modeled as a spring-damping system with a fixed saturation constraint to avoid the electro-mechanical generator nominal work area. This kind of control logic allowed us to define the PTO torque in Equation (17):

$$T_{PTO}=\left\{\begin{array}{cc}k\epsilon +c\dot{\epsilon }\hfill & \begin{array}{cc}& \mathrm{if}-T_s\le k\epsilon +c\dot{\epsilon }\le T_s\mathrm{and}-P_s\le k\epsilon \dot{\epsilon }+c{\dot{\epsilon }}^2\le P_s,\hfill \end{array}\hfill \\ T_s\hfill & \begin{array}{cc}& \mathrm{if}k\epsilon +c\dot{\epsilon }>T_s\mathrm{and}-\frac{P_s}{T_s}\le \dot{\epsilon }\le \frac{P_s}{T_s},\hfill \end{array}\hfill \\ -T_s\hfill & \begin{array}{cc}& \mathrm{if}k\epsilon +c\dot{\epsilon }<-T_s\mathrm{and}-\frac{P_s}{T_s}\le \dot{\epsilon }\le \frac{P_s}{T_s},\hfill \end{array}\hfill \\ \frac{P_s}{\dot{\epsilon }}\hfill & \begin{array}{cc}& \mathrm{if}k\epsilon \dot{\epsilon }+c{\dot{\epsilon }}^2>P_s\mathrm{and}\dot{\epsilon }\ge \frac{P_s}{T_s}\mathrm{or}\dot{\epsilon }\le -\frac{P_s}{T_s},\hfill \end{array}\hfill \\ -\frac{P_s}{\dot{\epsilon }}\hfill & \begin{array}{cc}& \mathrm{if}k\epsilon \dot{\epsilon }+c{\dot{\epsilon }}^2<-P_s\mathrm{and}\dot{\epsilon }\ge \frac{P_s}{T_s}\mathrm{or}\dot{\epsilon }\le -\frac{P_s}{T_s}.\hfill \end{array}\hfill \end{array}\right.$$

(17)

Parameters k and c represent stiffness and damping coefficients, respectively. $T_s$ corresponds to the saturation torque, and $P_s$ the saturation power. The ISWEC system employs a “slow-tuning” approach [54], adjusting control parameters such as the flywheel rotational velocity $\dot{\psi }$, c, and k in response to the spectral characteristics of the incoming sea state to maximize the energy harvested. This method simplifies control compared to wave-by-wave strategies, optimizing power production by adapting parameters for the entire duration of a sea state determined through its PSD. The “slow-tuning” concept signifies parameter adjustments based on evolving sea-state spectral properties while disregarding instantaneous wave profiles.

3.2. Framework

A classic Monte Carlo process was adopted with the aim of propagating uncertainties and avoiding assumptions inherent in the procedure. Two variables were defined as being subject to uncertainty. The first parameter was the device’s pitch inertia $M_{(5,5)}$. Previous analyses [45] highlight its impact on device performance. To propagate uncertainty concerning this variable, it was necessary to define a probability distribution for the associated uncertainties. A uniform distribution within the probability box defined by the interval [$-10\%M_{(5,5),0}$, $+10\%M_{(5,5),0}$] was chosen for the uncertainty of $M_{(5,5)}$, where $M_{(5,5),0}$ is the nominal expected value for the parameter. The choice was made after considering the effects on the pitch inertia of the device due to variations in the distribution of masses of the stern and bow ballasts. Consistent with the considerations for the first parameter described above, the second parameter subjected to real-world uncertainties was the position along the z-axis of the device’s center of gravity $Co{G}_z$. For this parameter, the uncertainty was modeled as an additive variation from the nominal value $Co{G}_{z,0}$: [$Co{G}_{z,0}-{\Delta }_{Co{G}_z}$,$Co{G}_{z,0}+{\Delta }_{Co{G}_z}$], describing again a uniform distributed probability-box. Different values of ${\Delta }_{Co{G}_z}$ additive noise were set, and their relative impact in the robustness quantification was observed. The output of the robustness analyses was the normalized power ($P_{out,adim}$) generated through the device via harvesting wave energy, denoted as $P_{out}$, normalized with respect to the maximum generated power in the whole dataset of the analyzed device $P_{out,adim}^{max}$:

$$P_{out,adim}=\frac{P_{out}}{P_{out}^{max}}.$$

(18)

While the literature includes studies on optimizing sampling in terms of quantity, distribution, and parameter space geometry [49], the focus for this study was solely calibrating the number of samples. To determine the minimum number of samples necessary to achieve a good system PDF representation, the focus was the convergence of standard deviation, as illustrated in Figure 4; that is, a good convergence ratio was reached with 1000 samples. The notation symbol $^$ represents the normalization of the nominal expected value of the output under examination, $f_m(\overrightarrow{x_i},\overrightarrow{u_0})$.

Simultaneously, when observing that the trend of the standard deviation converged with an increasing number of samples, a significant rise in computational costs could also be noticed. This is visually reported in Figure 5: the computational burden follows a likely linear trend, also highlighting the impracticality of investing such a substantial amount of time (8 h) to conduct a complete robustness quantification process for a single device. This drawback becomes even more relevant when considering the necessity of repeating a high number of simulations, e.g., for an RDO framework, for which the evaluation of the objective function might need to be performed thousands of times for the RQ phase. To mitigate this challenge and maintain a stochastic RQ approach, the next step was to develop a suitable surrogate model that was also able to reduce the computational costs for the uncertainty propagation process. A surrogate model, also known as a metamodel, can take the form of a mathematical relationship or an algorithm that captures the relationships between input and output variables; i.e., the procedure typically involves studying the output and input relationships and then fitting the right metamodels to represent that behavior. That is, a surrogate model $\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)$ approximates the computationally expensive model $\mathcal{M}\left(\overrightarrow{\upsilon }\right)$ [57].

Highlighting the crucial aspects of achieving an efficient metamodel, it becomes possible to list the characteristics of what we define as an ideal surrogate model (a gold model, as referred to in [57]). The surrogate model we obtained is superior not only if it is characterized by good accuracy but also if it has the ability to be trained with a limited number of computationally expensive model samples, i.e., saving computational time.

The technique chosen was Gaussian process regression (GPR), which can be considered an extension of the Kriging method [58] and was already used in [59], where a data-driven and GPR-based approach was developed in order to address the optimization problem for the PTO control strategy of the Pendulum Wave Energy Converter. Considering a dataset containing $N_{sample}$ joint observations of the extracted normalized power $y_p$ and uncertain parameters ${\overrightarrow{\upsilon }}_p\in {\mathbb{R}}^b$, the following equation was used:

$$\mathcal{D}=\{({\overrightarrow{\upsilon }}_p,y_p):p=1\dots N_{samples}\}.$$

(19)

In the present study, $b=2$ and each ${\overrightarrow{\upsilon }}_p$ represents a pair of $M_{(5,5)}$ and $Co{G}_z$. Our goal was to develop a metamodel $\tilde{\mathcal{M}}$, assuming that the observed extracted normalized power values followed a joint Gaussian distribution perturbed by an additional Gaussian noise $\zeta$:

$$\mathcal{M}\left(\overrightarrow{\upsilon }\right)=\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)+\zeta ,\zeta \sim \mathcal{N}(0,\lambda ).$$

(20)

Here, $\mathcal{M}\left(\overrightarrow{\upsilon }\right)$ is the computational expensive model underlying the observations $y_i$, and $\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)$ is a Gaussian process defined as follows:

$$\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)\sim GP(E\left(\overrightarrow{\upsilon }\right),\mathsf{\Omega }(\overrightarrow{\upsilon },\overrightarrow{{\upsilon }^{\prime }})).$$

(21)

This Gaussian process is characterized by the mean function $E\left(\overrightarrow{\upsilon }\right)$ and covariance function $\mathsf{\Omega }(\overrightarrow{\upsilon },\overrightarrow{{\upsilon }^{\prime }})$. The fundamental concept behind GPR is that any set of observations has a joint normal distribution with expected values provided via $E\left(\overrightarrow{{\upsilon }_p}\right)$ and a covariance matrix with entries represented by $\mathsf{\Omega }(\overrightarrow{{\upsilon }_p},\overrightarrow{{\upsilon }_j^{\prime }})$ for $p,j=1\dots N_{sample}$. In practical terms, it is possible to designate $f_p\left({\overrightarrow{\upsilon }}_p\right)=\tilde{\mathcal{M}}\left({\overrightarrow{\upsilon }}_p\right)$ as the observed values of the Gaussian process (distinct from $y_p$ due to noise subtraction) and ${\tilde{\mathcal{M}}}_{\ast }\left({\overrightarrow{\upsilon }}_{\ast }\right)$ as a random variable associated with the value of interest $f_{\ast }\left({\overrightarrow{\upsilon }}_{\ast }\right)$, signifying the prediction for ${\overrightarrow{\upsilon }}_{\ast }$. Therefore, it is possible to compute the model prediction as the mean value of ${\tilde{\mathcal{M}}}_{\ast }\left({\overrightarrow{\upsilon }}_{\ast }\right)$. In [60], a formal demonstration of the GPR mathematical aspects can be found.

To test the GPR approach’s accuracy and efficiency, a comparison with an expensive Monte Carlo simulation with 1000 samples was performed. Two different metamodels were tested with different expensive observations (50 and 100), and Latin Hypercube Sampling (LHS) was chosen as a sampling method. LHS enables the creation of an N-dimensional grid, with its density varying based on the desired number of samples. In the present case of study, N remained constant at 2. This grid was comprised of a set number of cells, corresponding to the total number of desired samples, ensuring that each cell was evaluated with a sample. This approach prevents the potential concentration of samples in specific regions of the parameter space (oversampling) or the under-sampling of certain areas, resulting in a uniform and well-distributed sampling pattern.

In Figure 6, the two obtained surfaces are represented, and in

Table 3. Comparison of an expensive Monte Carlo model $\mathcal{M}\left(\overrightarrow{\upsilon }\right)$ simulation and a GPR surrogate model $\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)$ simulation built with 50 and 100 samples, with an associated relative error ${\epsilon }_{rr,50}$${\epsilon }_{rr,100}$.

Index	$\mathcal{M}\left(\overrightarrow{\mathit{\upsilon }}\right)$	${\tilde{\mathcal{M}}}_{50}\left(\overrightarrow{\mathit{\upsilon }}\right)$	${\tilde{\mathcal{M}}}_{100}\left(\overrightarrow{\mathit{\upsilon }}\right)$	${\mathit{\epsilon }}_{\mathit{rr},50}$	${\mathit{\epsilon }}_{\mathit{rr},100}$
NRMSE [/]	—	$0.03$	$0.02$	—	—
Time [min]	$527.66$	$27.00$	$52.80$	—	—
$\widehat{\mu }$ [/]	$1.0811$	$1.0784$	$1.0821$	$0.25$ [%]	$0.09$ [%]
$\widehat{\sigma }$ [/]	$0.2760$	$0.2798$	$0.2743$	$1.38$ [%]	$0.62$ [%]
$sup\widehat{f}$ [/]	$0.6456$	$0.6216$	$0.6180$	$3.72$ [%]	$4.28$ [%]
$\widehat{\mu }+k\widehat{\sigma }$ [/]	$2.7369$	$2.7570$	$2.7279$	$0.73$ [%]	$0.33$ [%]
R [/]	$0.3687$	$0.3690$	$0.3684$	$0.08$ [%]	$0.08$ [%]
Q [/]	$0.7219$	$0.7230$	$0.7293$	$0.15$ [%]	$0.96$ [%]

the results are summarized. For the procedure, the surrogate model training and test datasets were different and without common samples. All the results were normalized with respect to the nominal expected value due to a nondisclosure agreement about the technology.

The outcomes show a clear overlap among the three surfaces. This means that the GPR potential is saturated; that is, increasing the number of observations resulted in a very similar surface. On the other hand, the accuracy of the metamodel can be considered good enough for statistical purposes. When looking at Figure 7, it is possible to reach the same conclusion. The distributions on the right (Figure 7b) differ for some localized values, but the general statistical behavior is similar among the two surrogates’ PDF and the expensive models. Moreover, this is consistent with the left image (Figure 7a) of the same Figure 7: the majority of data points align closely with the bisector of the plane defined by the real model’s normalized power output ($P_{out,adim}$—True on the y-axis) and the normalized power values predicted via surrogate models ($P_{out,adim}$—Pred on the x-axis). In the case of a perfectly accurate model, all data points would lie precisely on the bisector. Furthermore, there is evident minimal dispersion among the data points. The most significant deviations from the bisector (which describes the ideal prediction) are confined to a specific, limited region within the domain populated by only a few data points (<1%). The lack of samples made it difficult for the surrogate model to predict the output variable in those areas. When looking at Figure 7b, it can be seen that the green ($GP{R}_{100}$) and blue ($GP{R}_{50}$) curves do not follow the course of the red curve ($MonteCarl{o}_{1000}$) with the same accuracy as with the other intervals.

Another indicator of the good statistical approximation of the obtained surrogate model is provided via the value of the normalized root mean squared error (NRMSE). As suggested in [49], an NRMSE < $0.1$ is characterized by a reasonable predictive capacity. Thus, the surrogate model $\tilde{\mathcal{M}}\left(\overrightarrow{\upsilon }\right)$ is suitable for analysis, which would capture the overall system behavior, but accuracy in the prediction of specific points cannot be guaranteed. The trade-off between the speed and accuracy of the model is significantly directed towards the choice of ${\tilde{\mathcal{M}}}_{50}\left(\overrightarrow{\upsilon }\right)$ due to the metamodel accuracy, which is very similar to the one obtained via ${\tilde{\mathcal{M}}}_{100}\left(\overrightarrow{\upsilon }\right)$ with savings in computational costs equal to 25.8 min. The conservation of time resources is equal to 96% of the total amount of time requested for the RQ of a single device exploiting the ${\tilde{\mathcal{M}}}_{50}\left(\overrightarrow{\upsilon }\right)$ model: using 50 samples, it was possible perform the RQ for about twice the number of devices as would have been processed in the same time frame with ${\tilde{\mathcal{M}}}_{100}\left(\overrightarrow{\upsilon }\right)$ and with comparable and acceptable accuracy.

4. Results

The previously delineated framework was applied to the results of a nominal optimization process performed via a genetic algorithm; the selected site is near Pantelleria Island. In Figure 8, the scatter diagram for both the normalized annual sea energy density and occurrences is presented. Moreover, in the same Figure 8, the red dots represent the sea state set of representative waves (modeled as irregular waves according to a Jonswap spectrum) that were used in order to simulate the performances for each individual of the optimization algorithm. The waves were chosen to cover the relevant regions in terms of the energy and frequency of the scatter diagrams. The same representative sea states were employed in the RQ framework.

In Figure 9, the optimization results are reported in terms of $P_{out,adim}$ and $CapE{x}_{adim}$ (Equation (22)) as objective functions:

$$CapE{x}_{adim}=\frac{CapEx}{CapE{x}^{max}}.$$

(22)

The initial phase of this analysis centered on the examination of Pareto frontier devices with a direct comparison of outcomes following the assessment of previously expounded robustness indices for each WEC. Three different values of ${\Delta }_{Co{G}_z}$ were used to compare the responsiveness of the system robustness.

When examining Figure 10c, related to index $\widehat{\mu }$ in Figure 10, the results show that the majority of devices were located in the neighborhood of the normalized value of the mean performance ($\widehat{\mu }$) around 1, exhibiting good average behavior when subjected to the disturbances imposed on uncertain parameters.

As anticipated in the preceding sections, it was crucial to consider additional robustness metrics to obtain comprehensive information regarding the statistical trends of the examined performance.

When analyzing the $Q_i$ measurements in Figure 10b, it becomes evident that devices belonging to the nominal Pareto-set present PDF, in 90% of the occurrences, exceed values even lower than 65% of their relative $f_m(\overrightarrow{x_i},\overrightarrow{u_0})$, highlighting a difference between the claimed nominal performance and the performance achieved when taking into account uncertainties.

The achieved outcomes could be helpful during the decision-making phase. When a choice needs to be made in order to identify a singular device among the one belonging to the nominal Pareto set, the designer can rely on the robustness classification performed via the studied RIs.

The analysis now moves towards a broader pool of devices, which was selected from those involved in the whole nominal optimization process. The selection criteria were based on the normalized cost of energy ($Co{E}_{adim}$), defined in Equations (23) and (24) as the ratio between the device’s capital cost ($CapEx$ in $\left[MEuro\right]$) and the associated device’s annual energy production (AEP in $\left[\frac{MWh}{y}\right]$), considering an average device life time ($lif{e}_{25}$) of 25 years:

$$CoE=\frac{CapEx}{lif{e}_{25}AEP},$$

(23)

and normalized regarding the maximum dataset $CoE$:

$$Co{E}_{adim}=\frac{CoE}{Co{E}^{max}}.$$

(24)

The capital expenditures for the device $CapEx$ are composed of three primary factors: hull materials and construction, unit materials and construction, and PTO costs. Mooring system and operational and maintenance expenditures were excluded from the cost evaluation in this study due to their device-specific nature. The robustness of five thousand devices was quantified.

Also, in this case, the same three different ${\Delta }_{Co{G}_z}$ perturbations used above were applied. In the context of uncertainty propagation and result acquisition, the delineation of new Pareto frontiers occurred by establishing threshold values for each index. The limit values established above (for indices for which maximization indicated superior performance, Figure 11b,c,f) or below (for indices for which minimization indicated optimal performance, Figure 11a,d,e), for which robustness values were referred to, an individual WEC system’s performance can be considered acceptable. Therefore, these frontiers are not contingent upon the normalized nominal power output ($P_{out,adim}$) but are associated with the intrinsic robustness index pertinent to this parameter. From the results in Figure 11, it is evident that the imposed uncertainties play a significant role in shaping the Pareto set, considering the threshold set for the studied metrics. The most notable deviations were observed among the Pareto sets related to indices Q, $\widehat{\sigma }$, and $sup\widehat{f}$, respectively, in Figure 11b,e,f).

The standard deviation threshold was set to constrain its value equal to $5\%$ of the nominal expected value of the $i^{th}$ device. Moreover, the WCS for the singular system was bound to be equal to or greater than the $90\%$ of $f_m(\overrightarrow{x_i},\overrightarrow{u_0})$. The Q threshold was imposed to ensure that each perturbed WEC guaranteed a statistical trend with the $90\%$ of the occurrences defined by a $P_{out,adim}$ greater than the $90\%$ of the nominal expected values. In these three cases, the new robust-optimal devices exhibited significantly lower power levels compared to those belonging to the nominal Pareto set, revealing a trade-off between robustness and optimal performances. In Figure 11e, nearly 21% of the performance reduction could be achieved, and the $P_{out,adim}$ could decrease from 0.94–0.78 to 0.74–0.56 for Q and $sup\widehat{f}$, respectively, in Figure 11b,f. Moreover, imposing a threshold concerning device robustness led not only to a degradation in performance but also to a reduction in the Pareto set dimension itself, compelling the designer to make choices within a more limited range of objective functions ($P_{out,adim}$ and $CapE{x}_{adim}$). The mentioned Pareto size reduction is pronounced in Figure 11b,e,f.

Then, the study’s focus was directed toward emphasizing how the thresholds’ stringency magnitude impacts the robustness metric Q. Three distinct levels were investigated: 0.85 (the least stringent), 0.9 (medium-stringent), and 0.95 (the most rigorous). The outcomes are represented in Figure 12. Also, in this case, a trade-off between performance and robustness threshold severity was present.

At last, two different devices were selected from the whole pool of systems, namely D1 and D2, and a robustness performance assessment for the technology with respect to the site under analysis was carried out, i.e., Pantelleria Island (Figure 8). The case study was investigated in order to compare the performance robustness of the two different systems (D1 and D2), which both belong to the Pareto frontier and, therefore, are not dominant with respect to each other; i.e., neither of the two is univocally optimal compared to the other in relation to the performance examined ($CapE{x}_{adim}$, $P_{out,adim}$). The devices’ data are reported in

Table 4. D1 and D2 data.

Device	${\mathit{CapEx}}_{\mathit{adim}}$ [/]	${\mathit{P}}_{\mathit{out},\mathit{adim}}$ [/]	$\widehat{\mathit{\sigma }}$ [/]
D1	$0.45$	$0.74$	$0.13$
D2	$0.27$	$0.39$	$0.31$

.

The D1–D2 set of devices is characterized by a good trade-off between $CapE{x}_{adim}$ and performance for both systems: a medium-high normalized power for medium $CapEx$ for D1 and medium-low $P_{out,adim}$ but a low cost for D2. Although the two WECs are characterized by a large gap in terms of robustness, for those two devices, this value was quantified via the $\widehat{\sigma }$ index, i.e., $\Delta \sigma =0.18[/]$.

In Figure 13, the annual energy distribution over the normalized wave period and height is compared between D1 and D2. Here, light colors illustrate the variability in distribution when uncertainties were taken into account. On the x-axis, the wave period and significant height were normalized relative to the site’s most energetic period and wave height, denoted as the ratio between the relative sea state data and the site’s most energetic period ($T_{e,site}^{maxEn}$) and significant wave height ($H_{s,site}^{maxEn}$):

$$T_{e,adim}=\frac{T_e}{T_{e,site}^{maxEn}},$$

(25)

$$H_{s,adim}=\frac{H_s}{H_{s,site}^{maxEn}}.$$

(26)

The devices’ AEP was normalized relative to its respective distribution peak value ($AE{P}^{max}$):

$$AE{P}_{adim}=\frac{AEP}{AE{P}^{max}}.$$

(27)

The plot clearly shows the difference in terms of the standard deviation $\widehat{\sigma }$ for the two systems. The shading representative of the distribution fluctuation was wider for D2 than for D1. On the contrary, D1 showed a smaller range of variability. In this scenario, the decision-making process may favor the utilization of the D1 device due to its reduced susceptibility to uncertainties and higher productivity despite an increase in $CapE{x}_{adim}$ amounting to 18% of the maximum device cost in the whole dataset $CapE{x}^{max}$.

5. Conclusions and Future Works

The exploration of uncertain parameters’ impacts in real-world scenarios is in its early stages within the context of WEC design optimization processes. In the current literature addressing uncertainty treatment in WECs, the majority of attention is directed toward modeling control logics that are insensitive to uncertainties (i.e., robust control) and addressing uncertainties in the experimental phases of technology development.

When broadening the perspective to other research domains, it is observed that the most commonly employed methodologies in optimization processes taking into account perturbation are RO and RDO. Within the latter, particular emphasis is placed on robustness metrics aimed at quantifying the impact of parameter perturbations on the system’s performance.

The objective of this study was to investigate the information that various robustness indices can provide regarding optimal devices obtained following a nominal optimization process in the design of a WEC. For this purpose, a methodology based on stochastic uncertainty propagation via the Monte Carlo method was employed, exploiting a surrogate model (GPR) to reduce the computational cost of the RQ process.

The results highlight the necessity for attention to robustness evaluation in optimal WEC designs’ optimization process. To support this thesis, the results of the present study indicate that devices belonging to the nominal Pareto set are characterized by a $P_{out,adim}$ statistical trend that, in 90% of occurrences, was found to be higher than just 65% of the expected nominal value. Robustness indices were then evaluated for all individuals in a dataset obtained during the various generations of a design optimization process conducted using GA. These additional results reveal a significant deviation between the nominal Pareto set and the new one sought by setting differently stringent thresholds for robustness metrics. Around a 20% robustness performance reduction was achieved in the analyzed case when considering the standard deviation $\widehat{\sigma }$ index. Moreover, imposing a threshold on the device’s robustness can decrease the $P_{out,adim}$ from 0.94–0.78 to 0.74–0.56 for Q and $sup\widehat{f}$, respectively. These findings represent novel information that was not obtainable previously from the nominal optimization process, and it can constitute valuable data in decision-making. Furthermore, outcomes justified the intention to proceed toward the implementation of RO processes for WEC design.

Ultimately, a performance robustness assesment analysis was conducted for two selected devices chosen from the nominal Pareto set, namely D1 and D2. The aim of the case study was to delineate the annual productivity distribution over the sea $H_{s,adim}$ and $T_{e,adim}$. D1 and D2 were distinguished by a favorable balance between performance metrics and different degrees of robustness; i.e., the distribution fluctuation was wider for D2 than for D1. Therefore, the outcomes show that D1 WEC might be the better choice during a decision-making phase due to its reduced sensitivity to uncertainties and improved AEP despite a cost increase of around 18% of the maximum device $CapEx$ in the whole dataset.

The main limitations of this study were the high computational costs resulting from the necessary expensive numerical simulations to build the employed GPR model and the restriction of uncertain parameters to only two. The confidentiality constraints present in this study imposed restrictions concerning both the quantification of uncertainties’ impact on technology performance in absolute terms and the disclosure of devices’ design data. Nonetheless, the results obtained from analyzing different metrics effectively convey information to the designer about the robustness of different devices with respect to their expected performance ratings, and they allow for comparisons between different devices, providing insights into their PDF performance and robustness.

Moreover, the present research has revealed a certain trade-off between performance and robustness, which is not a new topic in the WEC literature. This sort of balance has been widely investigated in the robust control (e.g., WCS) and robust energy maximizing approaches [12,17]. Similarly, the performance compromise was investigated in WECs’ array field, e.g., with lower energy production costs, against the aggregate power performance degradation and the system sensitivity with respect to the wave interaction effects [61] or optimal array layout and efficiency [62].

Combining the obtained results and considering the limitations and trade-off highlighted in this study, future work, in addition to developing a WEC design RO framework, could involve a study of the impact that different types of control logic might have on the robustness of the conceptual design of a technology. Finally, the research should consider evaluating uncertainties in environmental and economic parameters, thus shifting attention towards the optimization of techno-economic indices such as the LCoE.