Optimal control of wave energy systems considering nonlinear Froude–Krylov effects: control-oriented modelling and moment-based control

Abstract

Motivated by the relevance of so-called nonlinear Froude–Krylov (FK) hydrodynamic effects in the accurate dynamical description of wave energy converters (WECs) under controlled conditions, and the apparent lack of a suitable control framework effectively capable of optimally harvesting ocean wave energy in such circumstances, we present, in this paper, an integrated framework to achieve such a control objective, by means of two main contributions. We first propose a data-based, control-oriented, modelling procedure, able to compute a suitable mathematical representation for nonlinear FK effects, fully compatible with state-of-the-art control procedures. Secondly, we propose a moment-based optimal control solution, capable of transcribing the energy-maximising optimal control problem for WECs subject to nonlinear FK effects, by incorporating the corresponding data-based FK model via moment-based theory, with real-time capabilities. We illustrate the application of the proposed framework, including energy absorption performance, by means of a comprehensive case study, comprising both the data-based modelling, and the optimal moment-based control of a heaving point absorber WEC subject to nonlinear FK forces.

1 Introduction

Wave energy converters (WECs) are devices designed to harvest the considerable energy available from ocean waves, by suitable conversion of the hydrodynamic energy contained in the surrounding wave field. Regardless of the specific conversion principle (see e.g. [12, 17, 56] for a detailed account of WEC absorption principles), it is already well-established that WEC systems intrinsically require suitable control technology to maximise energy absorption which, consequently, reduces the associated levelised cost of wave energy, hence directly supporting the pathway towards effective commercialisation [47, 67].

Among the desired features of WEC controllers, any candidate algorithm must respect three fundamental requirements, in order to be effective in realistic scenarios (see e.g. [18, 20]): (a) optimally maximise energy extraction from the wave resource, (b) minimise any potential risk of component damage (commonly formalised in terms of state and input constraints), and (c) perform the computation of the associated control law within the real-time limits characterising the specific WEC system. Since the fundamental purpose of a WEC is to effectively maximise energy absorption, which typically translates into an increase in the WEC operational range of motion [24, 47], nonlinear dynamical effects become significant, and non-representative linear models may become inaccurate [14]. While the computational effort usually increases considerably with the complexity associated with the nonlinearities included in the model, the corresponding gain in accuracy depends on the relevance of each nonlinear effect for the particular device under scrutiny. It is hence straightforward to see that a potential conflict arises between these three requirements, since achieving (a) and (b) above often require a precise (and hence potentially complex) model of the WEC dynamics/control objective, which almost inevitably conflicts with (c). In other words, a more accurate model/control objective representation leads to an increase in required computational (and analytical) effort to solve for the corresponding optimal control law which, in turn, can preclude real-time implementation/feasibility.

Though the vast majority of WEC control strategies, available in the literature, consider linear hydrodynamic WEC models (see e.g. [20]), some exceptions do exists, most of which consider relatively ‘simple’ (from an analytical complexity perspective) nonlinear hydrodynamic effects, including viscous drag forces [4, 7, 27], and so-called nonlinear (state-dependent) restoring effects [27, 49]. Nonetheless, for a large variety of devices currently in development, such as heaving point absorber WEC systems [36], a significant nonlinear contribution of the hydrodynamic force is the so-called Froude–Krylov (FK) effect (or force), which directly arises as the integration of the incident pressure field over the wetted surface of the device [33]. As a matter of fact, recent WEC design trends (see e.g. [11]) are based upon floating structures with a variable cross sectional area, and hence nonlinear FK effects can become potentially dominant. Though a significant effort has been expended in accurate numerical modelling of FK forces, as reported in the WEC literature (for both static and dynamic cases), e.g. [37, 39, 58], optimal control design and synthesis, capable of effectively taking into account such nonlinear FK effects, is significantly less prevalent, with any notable exceptions listed in the following paragraph.

The authors of [16] design a variable-structure (sliding mode, in this case, see e.g. [77]) control strategy for a heaving point absorber WEC, under both static, and dynamic FK forces, while [50] proposes a nonlinear model predictive controller for a WEC system with very similar dynamic characteristics. We note that both strategies presented in [16] and [50] share one fundamental disadvantage, which automatically precludes direct implementation of such controllers in realistic scenarios: The analytical model utilised to represent FK forces (particularly dynamic FK effects) intrinsically assumes a regular (monochromatic) free-surface elevation, i.e. composed of a single frequency component. Note that this is a rather strong assumption (since real waves are panchromatic), with the definition of the FK model, and hence the subsequent design and synthesis procedure for each specific control strategy, depending upon the availability of quantities such as the so-called wave number (see e.g. [29, 47]), which can only be defined for regular wave inputs, hence precluding a direct application to stochastic (irregular) wave fields. An analogous issue is present in the studies [36, 46], which implement so-called latching control of heaving WEC systems subject to nonlinear FK effects, by assuming that the free-surface elevation is effectively regular, so as to be able to compute a closed-form expression for the so-called latching time (which is ill-posed in the case of irregular wave inputs¹).

The lack of a model-based optimal control design for WEC systems, subject to nonlinear FK effects, can be attributed, at least partially, to the apparent unavailability of control-oriented models, representing such effects in a form ‘compatible’ with state-of-the-art control techniques. Though highly efficient numerical schemes have been presented in the literature of hydrodynamic WEC modelling, the assumptions adopted are often restrictive, and produce mathematical representations which are not entirely suitable for control design/synthesis. For instance, a particularly numerically efficient approach to the computation of nonlinear FK effects, is that proposed in [35, 36, 39], implemented via the open-source toolbox Nlfk4all [32, 34]. While such an approach is effectively able to achieve real-time computation of FK forces, the methodology proposed assumes a ‘frequency-by-frequency’ decomposition of the associated pressure field, an assumption which renders a mathematical description generally incompatible with model-based control design and synthesis procedures, which require a closed-form of (at least) the input-output description of the WEC dynamics [20].

Motivated by the lack of both suitable control-oriented modelling techniques for FK effects, and optimal control algorithms capable of effectively incorporating such nonlinear behaviour into the computation of a corresponding energy-maximising control law, we establish, in this paper, two main objectives. Firstly, we propose an exclusively data-based framework for the approximation of nonlinear FK effects in terms of mathematical structures compatible with state-of-the-art control procedures. We achieve this by computing a representative set of FK output data, via the numerical solver Nlfk4all, and subsequently utilising suitable techniques arising from the field of system identification and model reduction (see e.g. [1, 71]), tailored for each specific nonlinear FK effect, i.e. static and dynamic. Secondly, we propose a nonlinear optimal control strategy capable of effectively incorporating the computed data-based control-oriented description of FK forces into the energy-maximising control formulation, by exploiting the system-theoretic notion of a moment [2, 3, 28]. Moments are mathematical objects directly connected to the steady-state response of the WEC system, and have been recently shown to be highly efficient in parameterising the WEC optimal control problems (OCPs) in, for example, [23, 27]. In particular, [27] shows that the moment-based direct transcription of the WEC energy-maximising OCP renders a finite dimensional nonlinear program with appealing characteristics, i.e. uniqueness of the corresponding control parameterisation, and existence of globally optimal solutions, which facilitates the utilisation of efficient optimisation solvers with real-time capabilities. Given that the framework presented in [27] does not explicitly consider nonlinear FK effects, we present, in this paper, an extension of such a moment-based control methodology to the case of WEC systems containing FK nonlinearities, being able to preserve the attractive properties characterising the moment-based direct transcription of [27].

To briefly summarise, this paper encompasses the following three main contributions:

• (1) A data-based algorithm, tailored for control- oriented modelling of static nonlinear FK effects, based upon techniques belonging to the field of model reduction.

• (2) A data-based algorithm, tailored for control- oriented modelling of dynamic nonlinear FK effects, based upon techniques arising in the field of system identification.

• (3) A nonlinear optimal energy-maximising control algorithm based upon the system theoretic notion of moments, capable of effectively incorporating the models computed via (1) and (2), resulting in a well-posed optimisation problem, i.e. with guarantees of uniqueness of the proposed control parameterisation, and existence of globally optimal solutions.

We demonstrate, featuring an extensive case study, that the proposed data-based modelling procedure fits well with the control design requirements, i.e. there is a well-defined synergy between contributions (1), (2) and (3) of this paper, resulting in an overall control procedure capable of achieving maximum energy extraction from the wave resource, with real-time capabilities. This, in turn, represents a powerful tool for optimising operation of a wide range of WEC technology via tailored model-based control, including devices following recent design trends with variable cross-sectional area. Furthermore, we note that contributions (1) and (2) effectively constitute a systematic methodology to compute control-oriented WEC models incorporating nonlinear FK effects, and can also be used for by a general class of (alternative) model-based control formulations, such as nonlinear model predictive control (see [20]), further highlighting the potential impact of our study beyond the proposed moment-based controller (3).

The remainder of this paper is organised as follows. Section 1.1 introduces the key notations/conventions utilised throughout our paper, for clarity. Section 2 presents the fundamentals behind WEC hydrodynamic modelling, including explicit definitions for both static, and dynamic, nonlinear FK sources. Section 3 presents the complete proposed data-based control-oriented modelling procedure, while Sect. 4 discusses and derives the theoretical tools required to solve the optimal control problem for WEC systems with nonlinear FK effects using moment-based theory. Section 5 offers a case study which integrates both main contributions of this paper, i.e. data-based modelling and optimal control of WEC systems under nonlinear FK forces, for a heaving point absorber device, demonstrating the capabilities of the proposed framework. Finally, Sect. 6 encompasses the main conclusions of our study.

1.1 Notation/conventions

We consider standard notation throughout the reminder of this paper, with any exception detailed in the following. Note that we classify the considered notation according to its nature/functionality, for the sake of clarity.

1.1.1 Sets

$\mathbbm {R}^+$ ($\mathbbm {R}^{-}$) denotes the set of non-negative (non-positive) real numbers. $\mathbbm {C}^{0}$ denotes the set of pure-imaginary complex numbers, while $\mathbbm {C}_{<0}$ denotes the set of complex numbers with negative real-part. The notation ${\mathbb {N}}_{q}$ indicates the set of all positive natural numbers up to q, i.e. ${\mathbb {N}}_{q} = \{1,2,\ldots ,q\}$.

1.1.2 Scalars, vectors and matrices

The symbol 0 stands for any zero element, dimensioned according to the context. The symbol ${\mathbb {I}}_{n}$ denotes the identity matrix in $\mathbbm {C}^{n\times n}$, while the notation $\mathbf {1}_{n \times m}$ is used to denote a $n \times m$ Hadamard identity matrix (i.e. a $n \times m$ matrix with all its entries equal to 1). The superscript $\intercal$ denotes the transposition operator. The spectrum of a matrix $A\in {\mathbb {R}}^{n\times n}$, i.e. the set of its eigenvalues, is denoted by $\lambda (A)$. The symbol $\bigoplus$ denotes the direct sum of n matrices, i.e. $\bigoplus _{i=1}^{n} A_i = \text {diag}(A_1,A_2,\ldots ,A_n)$. The notation $\Re (z)$ and $\Im (z)$, with $z\in \mathbbm {C}$, stands for the real-part and the imaginary-part operators, respectively. The Kronecker product between two matrices $M_1\in \mathbbm {R}^{n \times m}$ and $M_2\in \mathbbm {R}^{p \times q}$ is denoted by $M_1\otimes M_2$ $\in \mathbbm {R}^{np \times mq}$. The vectorisation operator acting on a matrix $A\in \mathbbm {C}^{n\times m}$ is denoted as $\text {vec}(A) = [A]\in \mathbbm {C}^{nm}$. The symbol $\varepsilon _n\in \mathbbm {R}^{n}$ denotes a vector with odd entries equal to 1 and even entries equal to 0.

1.1.3 Functions

The generalised Dirac-$\delta$ function is denoted as $\delta :\mathbbm {R}\rightarrow \mathbbm {R}$, $t\mapsto \delta (t)$. Given two functions $f_1$ and $f_2$, such that $f_1:\mathscr {X} \rightarrow \mathscr {Y}$ and $f_{2}:\mathscr {Z} \rightarrow \mathscr {X}$, the composition $f_1(f_2(z))$, which maps all $z\in \mathscr {Z}$ to $f_1(f_2(z))\in \mathscr {Y}$, is denoted as $f_1 \circ f_2$. The convolution between two functions $g_1$ and $g_2$ over $\mathbbm {R}$, i.e. $\int _\mathbbm {R}g_1(\tau )g_2(t-\tau )\mathrm{d}\tau$ is denoted as . Let $w_1$ and $w_2$ be two functions belonging to the space $L^2(\varXi )$, with $\varXi \subset \mathbbm {R}$ closed. Then, the inner-product between $w_1$ and $w_2$ is given by $\langle w_1,w_2 \rangle = \int _{\varXi }w_1(\tau )w_2(\tau ) \mathrm{d}\tau$.

2 Hydrodynamic WEC modelling

In this section, we recall fundamental concepts behind nonlinear hydrodynamic modelling for wave energy conversion systems, based on potential flow theory (see e.g. [48, 55]). From now on, we assume a single degree-of-freedom (DoF) device, both for simplicity of exposition, and clarity of notation. However, note that we do this without any loss of generality, since the discussion on these theoretical preliminaries, the corresponding data-based approximation framework (proposed in Sect. 3), and moment-based control solution (derived in Sect. 4), can be extended to multi-DoF systems by following analogous procedures.

Let $z:\mathbbm {R}^+\rightarrow \mathbbm {R}$, $t\mapsto z(t)$ and $\eta : \mathbbm {R}^+ \rightarrow \mathbbm {R}$, $t\mapsto \eta (t)$, be the device excursion (displacement), and undisturbed free-surface elevation (measured at the centre of the body’s reference frame), respectively. The equation of motion of such a WEC system can be described in terms of a system $\varSigma _{\text {W}}$ written, for $t\in \mathbbm {R}^+$, [15, 29, 37] as

$$\begin{aligned} \varSigma _{\text {W}}: \left\{ \begin{aligned} m\ddot{z}&= f_{\text {r}}(z) + f_{\text {v}}({\dot{z}}) + f_{\text {d}}(\eta ) + \\&\quad \quad \, f^{\text {st}}_{\text {FK}}(\eta ,z) + f^{\text {dyn}}_{\text {FK}}(\eta ,z) - f_{\text {u}},\\ y&= C\begin{bmatrix} z&{\dot{z}} \end{bmatrix}^{\intercal }, \end{aligned} \right. \end{aligned}$$

(1)

where $m\in \mathbbm {R}^+$ is the generalised mass of the device, $f_{\text {r}}:\mathbbm {R}\rightarrow \mathbbm {R}$ is the radiation force, $f_{\text {v}}:\mathbbm {R}\rightarrow \mathbbm {R}$ represents viscous effects, $f_{\text {d}}:\mathbbm {R}\rightarrow \mathbbm {R}$ is the diffraction force, and the mappings $f^{\text {st}}_{\text {FK}}: \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$ and $f^{\text {dyn}}_{\text {FK}}: \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$, represent the so-called static and dynamic Froude–Krylov (FK) effects (or forces), respectively. Note that such FK forces, which constitute the main concern of this study, are subsequently described in detail in Sect. 2.1. The output $y:\mathbbm {R}^+\rightarrow \mathbbm {R}$ is assumed to be a linear combination of device displacement and velocity, defined via the matrix $C^{\intercal }\in \mathbbm {R}^2$.

Remark 1

Common (i.e. standard) choices for C in (1) include $C_{z} = [1\,\, 0]$ and $C_{{\dot{z}}} = [0\,\, 1]$, which set either displacement or velocity as the output y of system (1), respectively.

The radiation force $f_{\text {r}}$ is modelled based on linear potential theory and, using the well-known Cummins’ equation [13], can be written as

(2)

where $m_{\infty }$ is the so-called added-mass at infinite frequency, and $k_{\text {r}}:\mathbbm {R}^+ \rightarrow \mathbbm {R}$, $k_{r }\in L^2(\mathbbm {R})$, is the (causal) radiation impulse response function describing the memory effect of the fluid response.

Remark 2

The impulse response function $k_{\text {r}}$ fully characterises a linear time-invariant system with input ${\dot{z}}$ and output $f_{\text {r}}$. The rationality behind representing $f_{\text {r}}$ in terms of a convolution operator stems from the fact that the impulse response $k_{\text {r}}$ is virtually always computed numerically in a non-parametric form, via so-called boundary element methods (BEMs), see e.g. [5, 57]. We note that standard (in the sense of system dynamics) finite-dimensional parametric forms associated with $k_{\text {r}}$ can be computed via system identification procedures, see e.g. [22, 42, 59], though we do not pursue such an approximation in this paper, since the proposed control solution, presented in Sect. 4, can effectively handle (2) without the need of a specific parametric representation.

The mapping $f_{\text {v}}$, which represents viscous effects, is written in terms of a smooth approximation of the so-called Morison equation [53], i.e.

$$\begin{aligned} f_{\text {v}}({\dot{z}}) = -\alpha _{\text {v}} {\dot{z}} \sqrt{{\dot{z}}^2 + \varrho }, \end{aligned}$$

(3)

with $\varrho \in \mathbbm {R}^+$ sufficiently small, and $\alpha _{v}\in \mathbbm {R}^+$ directly depending on the geometric properties of the device.

The diffraction force, represented via $f_{\text {d}}$, can be described (analogously to the radiation force equation (2)), in terms of a convolution operator, i.e.

(4)

where the impulse response $k_{\text {d}}:\mathbbm {R}^+\rightarrow \mathbbm {R}$, $k_{\text {d}}\in L^2(\mathbbm {R})$, fully characterises a linear time-invariant system with input $\eta$ and output $f_{\text {d}}$. Note that, as in the case of the radiation force (see Remark 2), the diffraction kernel $k_{\text {d}}$ is virtually always computed numerically via BEM solvers, i.e. in a non-parametric form, hence the rational behind representing $f_{\text {d}}$ in terms of an appropriate convolution operator.

Finally, the map $f_{\text {u}}:\mathbbm {R}^+\rightarrow \mathbbm {R}$, $t\mapsto f_{\text {u}}(t)$, represents the control force (or law), which is to be designed so as to maximise the energy-absorption capabilities of the WEC system. The synthesis of such a control force, for the nonlinear WEC system defined in (1), effectively constitutes one of the main concerns of this study, and is specifically addressed, within a moment-based approach, in Sect. 4.

2.1 On the definition of nonlinear FK forces

As discussed in Sect. 1, a significant nonlinear component of the hydrodynamic force acting in (1) is the so-called FK effect (or force), which directly arises from integration of the incident pressure over the instantaneous wetted surface of the device. We revisit, in the reminder of this section, fundamental preliminaries associated with FK effects, as extensively discussed in, for instance, [36, 37, 39].

Let $\mathscr {S}_w(\eta , z) \equiv \mathscr {S}_w$ be the instantaneous wetted surface of the device, and let the mappings $p_{\text {st}}: \mathbbm {R}\rightarrow \mathbbm {R}$, and $p_{\text {dyn}}: \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$ denote the static and dynamic pressures, respectively, defined as

$$\begin{aligned} \begin{aligned} p_{\text {st}}(z) = \rho g z, \quad p_{\text {dyn}}(\eta , z) = \rho \frac{\partial \varPsi _{\text {I}}(\eta , z)}{\partial t}, \end{aligned} \end{aligned}$$

(5)

where $\varPsi _{\text {I}}$ denotes the so-called incident potential function (see [48, 55]).

The integration of the static pressure $p_{\text {st}}$, as defined in Eq. (5), over the instantaneous wetted surface $\mathscr {S}_w$, yields the so-called static FK force, i.e. the map $f^{\text {st}}_{\text {FK}}$, defined as

$$\begin{aligned} f^{\text {st}}_{\text {FK}}(\eta , z) = f_{\text {g}} - \iint _{\mathscr {S}_w} p_{\text {st}}(z) n_z \mathrm{d}\mathscr {S}, \end{aligned}$$

(6)

where $f_{\text {g}}$ denotes the gravity force, and $n_z$ is a (normalised) vector orthogonal to $\mathscr {S}_w$, while the corresponding integration of the dynamic pressure mapping $p_{\text {dyn}}$ in (5), i.e.

$$\begin{aligned} f^{\text {dyn}}_{\text {FK}}(\eta , z) = -\iint _{\mathscr {S}_w} p_{\text {dyn}}(\eta , z) n_z \mathrm{d}\mathscr {S}, \end{aligned}$$

(7)

gives rise to the so-called dynamic Froude–Krylov force $f^{\text {dyn}}_{\text {FK}}$.

In the case of a generic geometry, i.e. a general surface $\mathscr {S}_w$, the computation of both (6) and (7) can be performed by an appropriate spatial discretisation (e.g. mesh-based) of $\mathscr {S}_w$. Given the inherent complexity behind the numerical computation of the FK operators, defined in the paragraph above, analytical and semi-analytical solutions of both (6) and (7) have been proposed in [35, 36, 39]. In particular, the dynamic pressure mapping in (5) is explicitly written in terms of so-called Airy’s theory, while the instantaneous wetted surface is parameterised using different sets of coordinates, depending on the fundamental WEC geometry. As discussed in Sect. 1, while the approach proposed in [35, 36, 39], implemented in the open-source toolbox Nlfk4all [32, 34], is able to achieve real-time (numerical) computation of both static and dynamic FK effects, the methodology proposed assumes a ‘frequency-by-frequency’ decomposition of the associated pressure field, producing a mathematical description which is, in general, not compatible with state-of-the-art energy-maximising control design/synthesis procedures, whose formulation virtually always requires a closed-form of (at least) the input-output description of the WEC dynamics².

Motivated by the necessity of accurate, yet computationally efficient, control-oriented models of the WEC system in (1), we propose, in this paper, a data-based approximation framework for both static and dynamic FK forces to fulfil such a purpose. The approach, described in detail in Sect. 3, is capable of providing mathematical descriptions which are compatible with control design procedures, i.e. which are control-oriented, hence not only being convenient for the specific control approach presented in Sect. 4, but also for a general class of WEC control systems, such as those described in [20, 67].

3 Data-based control-oriented modelling of FK forces

Following the discussion provided in Sect. 2.1, we present, in this section, a control-oriented data-based approximation method for both static and dynamic FK mappings. To achieve such an objective, we adopt a system-theoretic approach to the definition of (6) and (7), and make explicit use of a specific type of deterministic persistently exciting (see e.g. [54]) signal termed a multisine (see e.g. [71, 72]), which we formally introduce in Sect. 3.1. Via the definition of such a class of signals, we collect representative output data for the systems induced by the mappings defined in (6) and (7), by means of the open-source nonlinear FK solver Nlfk4all. The corresponding input-output data set is then directly used in a black-box approach to identify suitable approximating structures via system identification procedures.

Given the rather different natures of the static and dynamic FK effects, the reminder of this section is organised into separate parts, as follows. Section 3.1 introduces the class of signals utilised for the approximation procedure proposed. Section 3.2 presents a general framework to approximate the static mapping of (6) in terms of the definition of a suitable function space. Section 3.3 discusses an approximation method for dynamic FK effects, in terms of a representative (see e.g. [14]) linear differential equation, i.e. a linear system designed to represent, as closely as possible, the associated FK mapping for the effective operational space of the device. Note that, as we explicitly discuss and demonstrate in Sect. 5, this is substantially different from linearising (7) about the device equilibrium position, which inherently derives a model designed to characterise the associated dynamic FK effects for an infinitesimal variation in displacement of the WEC device, which is at odds with WEC systems under optimal control conditions (see also the discussion provided in Sect. 1). Finally, Sect. 3.4 presents the overall control-oriented model for the WEC system under nonlinear FK effects, which is utilised for the optimal control design of Sect. 4.

3.1 Multisine excitation

A multisine signal is a specific type of deterministic excitation which can be used to solve a wide variety of identification problems [71, 72]. It has the advantage of being periodic, so that the issue of spectral leakage can be avoided. Furthermore, signals with user-controlled amplitude distribution and power spectrum can be easily produced, hence capable of exciting the target system in a pre-specified frequency band especially tailored for the wave energy modelling/control case³, due to the specific banded nature of the free-surface elevation (wave) input $\eta$. We give a formal definition of such signals in the following paragraph.

A multisine $f_{\text {id}}$ is a periodic signal with a band-limited spectrum. As such, it can be fully represented in terms of a Fourier series, i.e. a trigonometric sum of order K, with $K\in {\mathbb {N}}$ finite, such as

$$\begin{aligned} f_{\text {id}}(t) = \sum _{p=1}^K A_p \cos \left( l_p\omega _{\text {id}} t + \psi _p\right) , \end{aligned}$$

(8)

where $\{A_p, \psi _p\}_{i=1}^K\subset \mathbbm {R}$, and where $\omega _{\text {id}} = 2\pi /T_{\text {id}}$, with $T_{\text {id}}$ the measurement period, and $l_p$ a positive integer. The frequency $\omega _{\text {id}}$ is often referred to as the fundamental frequency (in [rad/s]) of the multisine signal. Typically, when a multisine signal is used, the set of phases $\{\psi _p\}$ is optimised aiming at minimising the so-called crest factor (see, for instance, [40]). One of the most well-established approaches, to achieve such an objective, is that proposed in [73], commonly referred to as Schroeder phases. In particular, for a flat amplitude spectrum (i.e. $A_p = C \in \mathbbm {R}^+$, $\forall p\in {\mathbb {N}}_K$), [73] suggests the set of phases in (8) as

$$\begin{aligned} \psi _p = \psi _1 - \dfrac{p(p-1)}{K}\pi , \end{aligned}$$

(9)

for $p\in {\mathbb {N}}_K/1$.

Remark 3

The definition of the set of Schroeder phases depends upon a rather simple analytical condition, i.e. equation (9). Although not considered in this manuscript, we do note that more ‘sophisticated’ methods exist for the optimal definition of the set of phases $\{\phi _p\}$, which are often based upon tailored optimisation routines. The reader is referred to, for instance, [40], for further discussion on this topic.

Fig. 1

Multisine excitation signal $f_{\text {id}}$: Time-domain representation (top), and associated frequency-domain amplitude spectra (bottom)

Aiming to briefly illustrate the nature of such a signal, Figure 1 (top) presents a multisine signal with measurement period $T_{\text {id}} \approx 314$ [s], which corresponds to a fundamental frequency of $\omega _{\text {id}} \approx 0.02$ [rad/s], together with an amplitude spectrum, as shown in Figure 1 (bottom). Note that the frequency band selected for the generation of such a multisine signal is set to [0.5, 5] [rad/s], which, as discussed in Sect. 5, is consistent with the identification procedures performed in the case study presented in this paper.

Remark 4

Note that this type of signal can excite a specific frequency band, with a user-defined spectrum (e.g. flat, as in the case of Figure 1), keeping an almost constant instantaneous amplitude in time. This last ‘quality’ is specifically useful for the wave energy case, facilitating an analogous definition of ‘regular wave height’, for the (inherently polychromatic) multisine signal $f_{\text {id}}$.

Remark 5

The limits associated with the exciting frequency band are intrinsically linked to the nature (i.e. dynamics), and operating conditions, of the system to be identified/approximated. We provide a more extensive discussion of the selection/design of such characteristics in Sect. 5.

3.2 Static FK effects

The mapping $f^{\text {st}}_{\text {FK}}$ is a static function which depends upon both the free-surface elevation $\eta$, and the displacement of the device z. Although, for some specific geometries, $f^{\text {st}}_{\text {FK}}$ can be derived analytically (see e.g. [36]), the objective of this section is to propose a ‘generic’ data-based methodology, i.e. independent of the specific device shape. To derive such a framework, we begin by noting that the mapping $f^{\text {st}}_{\text {FK}}$ can be seen, from a system-theoretic perspective, as a static system $\varSigma ^{\text {st}}$, with both $\eta$ and z as inputs, i.e.

$$\begin{aligned} \varSigma ^{\text {st}}: \left\{ y^{\text {st}} = f^{\text {st}}_{\text {FK}}(\eta , z),\right. \end{aligned}$$

(10)

where $y^{\text {st}}(t)\in \mathbbm {R}$ defines the corresponding output, i.e. the static FK force. It is relatively straightforward to see that such a system is inherently interconnected with the ‘remainder’ of the WEC dynamics described in (1). To be precise, one can define the auxiliary system $\varSigma _{W}/\varSigma ^{\text {st}}$ as

$$\begin{aligned} \varSigma _{\text {W}}/\varSigma ^{\text {st}}: \left\{ \begin{aligned} m\ddot{z}&= f_{\text {r}}(z) + f_{\text {v}}({\dot{z}}) + f_{\text {d}}(\eta ) + \\&\quad \, f^{\text {dyn}}_{\text {FK}}(\eta ,z) + f_{\text {u}},\\ y&= C_{z} \begin{bmatrix} z&{\dot{z}}\end{bmatrix}^{\intercal } = z, \end{aligned} \right. \end{aligned}$$

(11)

which is essentially system $\varSigma _{\text {W}}$ in (1) without considering static FK effects, and ‘decoupling’ $\varSigma ^{\text {st}}$ from (1), as illustrated in Fig. 2.

Remark 6

For the reminder of this section, and without any loss of generality, we consider $f_{\text {u}} = 0$, $\forall t\in \mathbbm {R}^+$. Note that $f_{\text {u}}$ does not play a role in the identification of (10), since $\varSigma ^{\text {st}}$ is fully characterised by a static mapping depending only directly on $\eta$ and z.

Fig. 2

Schematic illustration of static Froude–Krylov from a system-theoretic perspective

Before proceeding with a description of the corresponding data-based approximation framework for static FK effects, we introduce the following standing assumption: The mapping $f^{\text {st}}_{\text {FK}}$ belongs to the space generated by a family of real-valued functions $\{\phi _j\}_{j=1}^{\infty }$, with $\phi _j: \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$, $\phi _j\in {\mathcal {C}}^1$,$(\eta ,z)\mapsto \phi _j(\eta ,z)$, $\phi _j(0,0) = 0$, i.e. there exists a set of constants $\{a_j\}_{j=1}^{\infty }$ such that

$$\begin{aligned} f^{\text {st}}_{\text {FK}}(\eta , z) = \sum _{j=1}^{\infty }a_j \phi _j(\eta ,z), \end{aligned}$$

(12)

for every $\{\eta (t),z(t)\}\subset \mathbbm {R}$. This assumption, which is relatively standard in the literature (see the arguments posed in, for instance, [70]), provides a natural definition for an approximation of $f^{\text {st}}_{\text {FK}}$, as detailed in the following: From now on, we call the mapping ${\tilde{f}}^{\text {st}}_{\text {FK}}$, defined as

$$\begin{aligned} {\tilde{f}}^{\text {st}}_{\text {FK}}(\eta , z) = \sum _{j=1}^{N} a_j \phi _j(\eta ,z), \end{aligned}$$

(13)

with $N\in {\mathbb {N}}$ finite, the approximated static nonlinear FK force. This definition, posed by means of (13), is based upon the idea of ‘truncating’ the expansion for $f^{\text {st}}_{\text {FK}}$ in Eq. (12), considering only N basis functions, i.e. the static Froude–Krylov mapping is essentially approximated by its expansion on the subset $\mathscr {P} = \{\phi _j\}_{j=1}^N$. The main idea is now to compute the approximated FK mapping, in the sense of (13), solely based upon data, i.e. without explicit knowledge of the internal structure of (10).

Remark 7

For the sake of completeness, we note that, by means of (13), one can straightforwardly define an approximating system ${\tilde{\varSigma }}^{\text {st}}$, analogously to Eq. (10), simply as

$$\begin{aligned} {\tilde{\varSigma }}^{\text {st}}: \left\{ {\tilde{y}}^{\text {st}} = {\tilde{f}}^{\text {st}}_{\text {FK}}(\eta , z),\right. \end{aligned}$$

(14)

where, clearly, ${\tilde{y}}^{\text {st}} \approx y^{\text {st}}$ in the sense of (13).

Remark 8

In practice, the set of functions $\mathscr {P}$ can be selected via a trial and error procedure, starting with, for instance, a polynomial expansion, if the mapping $f^{\text {st}}_{\text {FK}}$ is known to be smooth.

In order to proceed with the proposition of the data-based approximation framework for static FK forces, we introduce the following set of auxiliary variables:

$$\begin{aligned} \begin{aligned}&{\mathcal {P}} = \begin{bmatrix} a_1&a_2&\ldots&a_N \end{bmatrix},\\&\varPhi (\eta , z) = \begin{bmatrix} \phi _1(\eta ,z)&\phi _2(\eta ,z)&\ldots&\phi _N(\eta ,z) \end{bmatrix}^{\intercal }, \end{aligned} \end{aligned}$$

(15)

where $\{{\mathcal {P}}^{\intercal }, \varPhi (\eta ,z)\}\subset \mathbbm {R}^N$, so that the approximated mapping defined in (13) can be written in a compact form as

$$\begin{aligned} {\tilde{f}}^{\text {st}}_{\text {FK}}(\eta , z) = {\mathcal {P}}\varPhi (\eta ,z). \end{aligned}$$

(16)

With the compact definition presented in (16), the approximation problem reduces to finding a suitable matrix ${\mathcal {P}}$, for any given basis-function vector $\varPhi (\eta ,z)$. In particular, we now propose a method to compute ${\mathcal {P}}$ based upon a recursive least-squares approach, inspired by the data-driven model reduction framework of [70]. Let ${\mathcal {T}}^w_{k} = \{t_{k-w-1},\ldots ,t_{k-1}, t_k\} \subset \mathbbm {R}^+$ be a set of time-instants, where we numerically evaluate the output of the target static FK system, i.e. $y^{\text {st}}$ in (10). Note that the set ${\mathcal {T}}^w_k$ basically represents a moving window of $w\in {\mathbb {N}}$ samples. Suppose that, in the numerical evaluation of $y^{\text {st}}$, we select the external input $\eta$ (i.e. free-surface elevation) as a band-limited multisine signal, i.e. $\eta = f_{\text {id}}$ (see Sect. 3.1), and let $\varXi _k\in \mathbbm {R}^{w\times N}$ and $\varUpsilon _k\in \mathbbm {R}^{w}$, with $w \ge N$, be defined as

$$\begin{aligned} \begin{aligned} \varXi _k&= \begin{bmatrix} \varPhi _{t_{k-w-1}}&\ldots&\varPhi _{t_{k-1}}&\varPhi _{t_{k}} \end{bmatrix}^{\intercal }, \\ \varUpsilon _k&= \begin{bmatrix} y^{\text {st}}_{t_{k-w-1}}&\ldots&y^{\text {st}}_{t_{k-1}}&y^{\text {st}}_{t_{k}} \end{bmatrix}^{\intercal }. \end{aligned} \end{aligned}$$

(17)

with $\varPhi _{t_r} = \varPhi (\eta (t_{r}), z(t_{r}))\in \mathbbm {R}^{N}$ and $y^{\text {st}}_{t_{r}} = y^{\text {st}}(t_{r})\in \mathbbm {R}$.

Let ${\mathcal {P}}_k$ be an online estimate of the matrix ${\mathcal {P}}$ in (15), obtained using the set $T^w_k$, i.e. computed at the time $t_k$ using the last w instants of time. The underpinning idea is to determine ${\mathcal {P}}_k$ based upon the following optimal criterion:

$$\begin{aligned} \min _{{\mathcal {P}}_k} \Vert \varXi _k{\mathcal {P}}_k^{\intercal } - \varUpsilon _k \Vert _2^{2}, \end{aligned}$$

(18)

i.e. to find the least-squares solution of $\varXi _k{\mathcal {P}}_k^{\intercal } = \varUpsilon _k$, in a recursive fashion. From now on, and aiming to guarantee well-posedness of the upcoming methodology, we assume that the elements of $T^w_k$ are selected such that $\varXi _k$ is full column rank for all $k\in {\mathbb {N}}$.

Remark 9

Note that, since $\eta$ is selected as a multisine signal (which is effectively a persistently exciting signal, see [54, 71]), the set $T^{w}_k$ can always be selected such that the above assumption holds, i.e. we do not lose any generality.

We now propose, in Algorithm 1, a recursive least-squares solution to solve for ${\mathcal {P}}_k$, by adapting the procedure proposed in [70], originally developed for model reduction purposes.

Algorithm 1 can be briefly summarised as follows. After selecting a suitable multisine signal $f_{\text {id}}$ for the generation of $\eta$, and a sufficiently large initial value $k^0$, the matrices $\varXi _k$ and $\varUpsilon _k$ are constructed iteratively, using information on the (supplied) input $\eta = f_{\text {id}}$, and the numerically computed device motion z, and total static force $y^{\text {st}}$. Note that the latter two variables can be readily computed via any numerical nonlinear FK solver, such as the Nlfk4all toolbox [32, 34] (which is explicitly considered in this paper within the case study presented in Sect. 5). Every time a new set of samples is available, old information is discarded, and the algorithm is repeated until achieving a certain threshold condition on the error between iterations, specified by the (sufficiently small) value $\epsilon$.

Algorithm 1 (Static FK approximation) Let $k^0$ be a sufficiently large integer, and let $\epsilon$ be a sufficiently small user-defined error tolerance. Define ${\mathcal {P}}^0\in \mathbbm {R}^{1\times N}$ as the initialisation vector for the computation of ${\mathcal {P}}$.

Remark 10

The matrix $\varPsi _k$ is always well-defined due to the fact that $\varXi _k$ is full column rank for all $k\in {\mathbb {N}}$.

Remark 11

The selection of ${\mathcal {P}}^0$, required to start the recursion, can be done in terms of ‘dummy’ values (e.g. random values with a uniform distribution). Since the associated recursive algorithm ‘forgets’ old measurements (i.e. the values involved in the computation of ${\mathcal {P}}$ are updated at each k), after a sufficient number of iterations the effect of such a selection is effectively forgotten.

3.3 Dynamic FK effects

In contrast to the case of static effects, the mapping associated with dynamic FK forces, i.e. equation (7), inherently represents a dynamic entity. In particular, taking a system-theoretic perspective, $f^{\text {dyn}}_{\text {FK}}$ can be generally represented as the output mapping $y^{\text {dyn}}$ of a dynamical system $\varSigma ^{\text {dyn}}$, with both $\eta$ and z as inputs, i.e.

$$\begin{aligned} \varSigma ^{\text {dyn}}: \left\{ \begin{aligned} {\dot{q}}&= \kappa (q,\eta ,z),\\ f^{\text {dyn}}_{\text {FK}}&= \theta (q,\eta ,z) = y^{\text {dyn}}, \end{aligned}\right. \end{aligned}$$

(19)

where $q(t)\in \mathbbm {R}^n$ is the associated state-vector, with $n\in {\mathbb {N}}$ sufficiently large. From a physical standpoint (see e.g. [33]), the (locally smooth) mappings $\kappa :\mathbbm {R}^n\times \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}^n$ and $\theta :\mathbbm {R}^n\times \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$ are such that the following conditions hold:

• 1. $\kappa (0,0,0) = 0$ and $\theta (0,0,0) = 0$.

• 2. The zero-equilibrium of system ${\dot{q}} = \kappa (q,0,0)$ is locally exponentially stable.

• 3. For every $(q,\eta ,z)$, defined in a sufficiently small neighbourhood ${\mathcal {N}}$ of $(q,\eta ,z) = (0,0,0)$, the behaviour of system (19) can be well-represented in terms of a strictly proper linear system ${\tilde{\varSigma }}^{\text {dyn}}$ with only $\eta$ as external input. To be precise, the following conditions

  $$\begin{aligned} \begin{aligned}&\left. \dfrac{\partial \kappa }{\partial x}\right| _{(0,0,0)} = F, \, \left. \dfrac{\partial \kappa }{\partial \eta }\right| _{(0,0,0)} = G, \, \left. \dfrac{\partial \kappa }{\partial z}\right| _{(0,0,0)} = 0, \\&\left. \dfrac{\partial \theta }{\partial x}\right| _{(0,0,0)} = H, \, \left. \dfrac{\partial \theta }{\partial \eta }\right| _{(0,0,0)} = 0, \,\,\, \left. \dfrac{\partial \theta }{\partial z}\right| _{(0,0,0)} = 0, \end{aligned}\nonumber \\ \end{aligned}$$

  (20)

  hold, with $F\in \mathbbm {R}^{n\times n}$, $\{G,H^{\intercal }\}\subset \mathbbm {R}^n$, and hence system ${\tilde{\varSigma }}^{\text {dyn}}$ can be written as

  $$\begin{aligned} {\tilde{\varSigma }}^{\text {dyn}}: \left\{ \begin{array}{l} {\dot{q}}(t) = Fq(t) + G\eta (t),\\ {\tilde{f}}^{\text {dyn}}_{\text {FK}}(t) = Hq(t + t_c) = {\tilde{y}}^{\text {dyn}}(t), \end{array}\right. \end{aligned}$$

  (21)

  where $t_c\in \mathbbm {R}^+$, and $\lambda (F)\subset \mathbbm {C}_{<0}$ by condition 2 above.

Remark 12

The value $t_c$, which denotes a non-causal time-shift (i.e. advance), stems from the fact that the generated wave (free-surface elevation) may impact the WEC body and exert a wave force before any wave has reached the device ‘centre’ (see, for instance, [29]).

There is, naturally, a significant motivation to work with the linear structure posed in (21). Such a representation is computationally simpler, obeys the principle of superposition, and lends itself to a vast set of mathematical tools which can be used for analysis, simulation, and control/estimator design. However, in the current literature, system ${\tilde{\varSigma }}^{\text {dyn}}$ is almost exclusively characterised via so-called boundary element method solvers, such as the open-source software Nemoh [5], where the impulse response associated with system (21) is computed (in a non-parametric form) under the assumption of infinitesimally small motion of the device about the zero-equilibrium (see [14, 18, 20]). Given that the design objective for WECs is that of maximising converted energy, which typically implies a large induced device motion (see also the discussion provided in Sect. 1), such a methodology is likely to result in a non-representative linear model for the dynamic FK effect.

In contrast to BEM solvers, we propose a methodology to compute a representative linear model ${\tilde{\varSigma }}^{\text {dyn}}$, valid for a given set of wave operating conditions for the device, i.e. significant wave heights and peak periods. To achieve this, we employ tools from the field of system identification, and we propose a framework to provide representative models via so-called black-box structures, using only input-output data in the frequency-domain. Such a methodology is discussed in the following paragraphs.

Let $\mathscr {U} = \{\eta _i = f_{\text {id}_i}\}_{i=1}^Q$ be a set of suitably selected multisine input signals (i.e. free-surface elevation profiles), with $Q\in {\mathbb {N}}_{\ge 1}$, generating a corresponding set of outputs $\mathscr {Y} = \{y_i^{\text {dyn}}\}_{i=1}^Q$ (i.e. dynamic FK forces). Let $(\eta _i, y_i^{\text {dyn}})$ denote an input-output pair of signals for system (21). We define the so-called empirical transfer function estimate (ETFE) $H_i:\mathbbm {C}^0\rightarrow \mathbbm {C}$, $j\omega \mapsto H_i(j\omega )$, for each input-output pair, in terms of the expression

$$\begin{aligned} H_i(j\omega ) = \dfrac{{\mathcal {N}}_i(j\omega )}{Y_i^{\text {dyn}}(j\omega )}, \end{aligned}$$

(22)

with $i\in {\mathbb {N}}_Q$, and where ${\mathcal {N}}_i:\mathbbm {C}^0 \rightarrow \mathbbm {C}$ and $Y_i^{\text {dyn}}:\mathbbm {C}^0 \rightarrow \mathbbm {C}$ denote the Fourier transform of $\eta _i$ and $y_i^{\text {dyn}}$, respectively. Via the set of ETFE mappings $\mathscr {H} = \{H_i\}_{i=1}^Q$, one can readily obtain the so-called average ETFE, ${\bar{H}}$, computed with the aim of building a low-variance set to use as input to the frequency-domain identification algorithm (see e.g. [63]). The explicit computation of ${\bar{H}}$ is simply given by:

$$\begin{aligned} {\bar{H}}(j\omega ) = \dfrac{1}{Q} \sum _{i=1}^Q H_i(j\omega ). \end{aligned}$$

(23)

Recall that the ultimate objective of the proposed system identification procedure is to obtain a parametric form which approximates the behaviour of $\varSigma ^{\text {dyn}}$ in terms of a representative linear structure ${\tilde{\varSigma }}^{\text {dyn}}$, as in equation (21), based upon the characterisation provided by the average ETFE (23), which is computed exclusively in terms of input-output data. Before discussing the method from an algorithmic perspective, and without any loss of generality, let us re-write the average ETFE in (23) as

$$\begin{aligned} {\bar{H}}(j\omega ) = e^{j\omega t_c}{\bar{H}}_c(j\omega ), \end{aligned}$$

(24)

where the term $e^{j\omega t_c}$ denotes the frequency-domain equivalent of the time-shift (advance) corresponding with the dynamic FK system (see Remark 12) and, hence, ${\bar{H}}_c$ only represents the causal component of ${\bar{H}}$.

The strategy considered here, to compute a state-space structure (as in (21)), from the average ETFE (23), is that of subspace-based identification [79]. In particular, we consider the methodology outlined in [78], which computes the associated Hankel matrices directly from frequency-domain data (i.e. the input-output information encoded in the mapping ${\bar{H}}$)⁴. Furthermore, we combine such a methodology with an iterative procedure to compute an estimate of the corresponding time-advance $t_c$. The proposed methodology is summarised in Algorithm 2, from a systematic perspective.

Algorithm 2 (Dynamic FK approximation) Let the input set $\mathscr {U} = \{\eta _i\}_{i=1}^Q$ be such that each $\eta _i$ is a suitably selected multisine signal $f_{\text {id}_i}$. Let ${\mathcal {T}}_c = \{t_{c_i}\}_{i=1}^P\subset [a_l, a_h]\subset \mathbbm {R}^+$ be a set of trial time-shifts, with $a_l$ sufficiently small, and $a_h$ and $P\in {\mathbb {N}}$ sufficiently large. Let ${\tilde{\varSigma }}$ denote an approximated state-space system, computed from an average ETFE ${\bar{H}}$, with a (finite) user-selected order (dimension) ${\tilde{n}}$, and let $\tilde{{\mathcal {H}}}$ denote its associated transfer function.

The strategy in Algorithm 2 can be briefly described as follows. Firstly, the user selects a suitable set of multisine input signals (see Sect. 3.1) and collects the corresponding dynamic FK output data for each defined free-surface elevation. The resulting input-output set is then used to construct the associated average ETFE as in (23). Secondly, and since the value for $t_c$ in (21) is, in general, unknown, a finite-set of trial values ${\mathcal {T}}_c$ is chosen. For each trial value contained in the set, ${\bar{H}}_c$ is constructed (as in equation (24)), and a corresponding approximating state-space structure ${\tilde{\varSigma }}$ of order (dimension) ${\tilde{n}}$ is computed using subspace techniques. Note that the latter step is indicated in Algorithm 2 with the ‘function’ identify$(\cdot \,, \cdot )$. Finally, the system ${\tilde{\varSigma }}_i$ producing the lowest fitting error ${\mathcal {E}}_i$, together with its associated $t_c$ value, is selected as the linear representative approximating model ${\tilde{\varSigma }}^{\text {dyn}}$, as in (21).

To finalise this section, and provide a corresponding overview of the approximation framework, we provide a schematic illustration of the final structure for the complete approximating system for nonlinear FK forces in Figure 3, including both static and dynamic effects. Note that, as depicted in Figure 3, the ultimate objective is to provide a (computationally and representationally) convenient approximation of the total FK force $y_{\text {FK}} = y^{\text {st}} + y^{\text {dyn}}$, via ${\tilde{y}}_{\text {FK}} = {\tilde{y}}^{\text {st}} + {\tilde{y}}^{\text {dyn}}$.

Fig. 3

Schematic representation of the target, and approximated, FK systems

3.4 Control-oriented model

With the proposed approximating structures for both static and dynamic FK effects, computed in Sects. 3.2 and 3.3, respectively, one can provide a control-oriented model ${\tilde{\varSigma }}_{\text {W}}$, approximating the nonlinear WEC dynamics $\varSigma _{\text {W}}$ in (1). To achieve such an objective, we take the steps detailed in the following paragraphs.

We begin by noting that, without any loss of generality, the mapping ${\tilde{f}}^{\text {st}}_{\text {FK}}$, characterising nonlinear static FK forces within the presented approximation framework, can be ‘separated’ as

$$\begin{aligned} {\tilde{f}}^{\text {st}}_{\text {FK}} = {\tilde{f}}^{\text {st}}_{\text {FK}_\text {l}} + {\tilde{f}}^{\text {st}}_{\text {FK}_{\text {nl}}}, \end{aligned}$$

(25)

where

$$\begin{aligned} \begin{aligned}&{\tilde{f}}^{\text {st}}_{\text {FK}_\text {l}}(z) = {\mathcal {P}}\left. \frac{\partial \varPhi }{\partial z}\right| _{(0,0)} z,\\&{\tilde{f}}^{\text {st}}_{\text {FK}_\text {nl}}(\eta ,z) = {\tilde{f}}^{\text {st}}_{\text {FK}}(\eta ,z) - {\tilde{f}}^{\text {st}}_{\text {FK}_\text {l}}(z), \end{aligned} \end{aligned}$$

(26)

i.e. into linear ${\tilde{f}}^{\text {st}}_{\text {FK}_\text {l}}$, and nonlinear ${\tilde{f}}^{\text {st}}_{\text {FK}_\text {nl}}$, contributions, respectively.

Remark 13

The linear map ${\tilde{f}}^{\text {st}}_{\text {FK}_\text {l}}$ depends only on z: From the physical laws governing the WEC dynamics, it is possible to show that $\left. \tfrac{\partial \varPhi }{\partial \eta }\right| _{(0,0)} = 0$, and hence the linear part of the static FK force depends only upon the device displacement. In fact, the coefficient ${\mathcal {P}}\left. \tfrac{\partial \varPhi }{\partial z}\right| _{(0,0)}$ acts exactly as the well-known hydrostatic stiffness coefficient, commonly used within linear potential flow theory (see e.g. [29]).

With the explicit nonlinear contribution of the approximated static FK force, made available via equations (25) and (26), it is possible to ‘condense’ the nonlinear characteristics of the WEC system (1) in terms of a ${\mathcal {C}}^1$ function $g_{\text {W}}:\mathbbm {R}\times \mathbbm {R}\times \mathbbm {R}\rightarrow \mathbbm {R}$, defined as

$$\begin{aligned} g_{\text {W}}(\eta ,z,{\dot{z}}) = f_{\text {v}}({\dot{z}}) + {\tilde{f}}^{\text {st}}_{\text {FK}_{\text {nl}}}(\eta ,z), \end{aligned}$$

(27)

which also includes the viscous effects, modelled via Eq. (3).

With respect to dynamic FK forces, and aiming to be consistent with the mathematical description considered for the diffraction force in (4), we write the output of the approximating system ${\tilde{\varSigma }}^{\text {dyn}}$ in terms of its associated impulse response function, i.e.

(28)

where ${\tilde{k}}^{\text {dyn}}_{\text {FK}}:\mathbbm {R}^+\rightarrow \mathbbm {R}$, ${\tilde{k}}^{\text {dyn}}_{\text {FK}}\in L^2(\mathbbm {R})$, can be directly defined in terms of the triple of matrices (F, G, H) composing (21) as

$$\begin{aligned} {\tilde{k}}^{\text {dyn}}_{\text {FK}}(t) = He^{F(t+t_c)}G. \end{aligned}$$

(29)

Using the results presented up until this point, the approximating WEC dynamical equation can be conveniently written as

(30)

with ${\mathcal {M}} = m + m_{\infty }$. Note that the linear and nonlinear contributions are explicitly separated in equation (30).

Remark 14

The function $k_{\text {d}} + k^{\text {dyn}}_{\text {FK}}$ resembles the so-called excitation force kernel, commonly employed within linear potential flow theory models, though with a fundamental difference: $k^{\text {dyn}}_{\text {FK}}$ is obtained here in terms of a representative linear model (see Sect. 3.3), as opposed to linear potential flow theory, where the corresponding model is computed under the assumption of infinitesimally small motion of the device about the zero-equilibrium .

Finally, and to provide a model compatible with control design (as approached in Sect. 4), let $x = [z,\,{\dot{z}}]^{\intercal }\in \mathbbm {R}^2$ be the state vector associated with the WEC system. The approximating control-oriented model can be then expressed as

(31)

with the map $g:\mathbbm {R}\times \mathbbm {R}^2\rightarrow \mathbbm {R}^2$ given by

$$\begin{aligned} g(\eta ,x) = B\,g_{\text {W}}(\eta , C_{z}x, C_{{\dot{z}}}x), \end{aligned}$$

(32)

and where the pair of matrices (A, B), defining (31), can be straightforwardly defined as

$$\begin{aligned} A = \begin{bmatrix} 0 &{} 1\\ {\mathcal {M}}^{-1}{\mathcal {P}}\left. \tfrac{\partial \varPhi }{\partial z}\right| _{(0,0)} &{} 0 \end{bmatrix},\quad B = \begin{bmatrix} 0 \\ {\mathcal {M}}^{-1} \end{bmatrix}. \end{aligned}$$

(33)

Note that the output of (31) is now fixed to be the device velocity, given its relevance within the definition of the associated energy-maximising optimal control problem for WECs (see Sect. 4 for further detail).

Remark 15

It is straightforward to show that $(x,\eta ,f_{\text {u}}) = (0,0,0)$ is an equilibrium point for (31). Furthermore, given that $\tfrac{\partial g(x,\eta )}{\partial x}\vert _{(0,0)} = 0$, with g as in (32), the local stability properties of the zero-equilibrium of ${\tilde{\varSigma }}_{\text {W}}$ are governed by the linearised equation

(34)

As it is well-known in the marine engineering community, the Volterra Eq. (34) can be shown to be asymptotically stable using dissipativity arguments, for any meaningful values of the parameters involved, and radiation kernel $k_{\text {r}}$ (see e.g. [29, 75]).

4 Moment-based optimal control

Based upon the model derived via the framework presented throughout Sect. 3, we now consider the design of a nonlinear moment-based optimal controller for a WEC system, subject to nonlinear FK effects. In particular, we present, in this paper, a generalisation of the moment-based control results presented in [27], by incorporating a general class of nonlinear wave-dependent effects within the computation of the corresponding optimal control law. Before getting into the specific details characterising the proposed control strategy, we formally introduce the energy-maximising OCP associated with the WEC control problem.

As discussed throughout Sect. 1, the control objective for wave energy systems is that of maximising energy extraction from the incoming wave field, via a suitably designed optimal control law $f^{\text {opt}}_{\text {u}}$, injected into the system via the so-called power take-off (PTO) system (i.e. actuator). To be precise, the associated OCP can be written in terms of an objective function ${\mathcal {J}}:\mathscr {U}_{\text {u}}\rightarrow \mathbbm {R}$, $u\mapsto {\mathcal {J}}(u)$, defined (see e.g. [18, 67]) as

$$\begin{aligned} {\mathcal {J}}(f_{\text {u}})= & {} \dfrac{1}{T}\int _{\varXi } f_{\text {u}}(\tau ){\dot{z}}(\tau ) \mathrm{d}\tau \nonumber \\ {}= & {} \dfrac{1}{T}\int _{\varXi } f_{\text {u}}(\tau )y(\tau ) \mathrm{d}\tau , \end{aligned}$$

(35)

where y is as defined in (31), $\varXi = [0,\,T]\subset \mathbbm {R}^{+}$, and $\mathscr {U}_{\text {u}}$ denotes the set of admissible inputs.

In addition to the control objective function posed in (35), the WEC OCP considers both state and input constraints, aiming to secure energy-maximisation from incoming waves, while guaranteeing the intrinsic safety limits associated with both WEC system (state), and actuator (input), components. Such a constraint set $\mathscr {C}$ can be formally written as

$$\begin{aligned} \mathscr {C}:\left\{ \begin{aligned}&\vert z(t) \vert \le Z_{\text {max}},\\&\vert {\dot{z}}(t) \vert \le V_{\text {max}},\\&\vert u(t) \vert \le U_{\text {max}}, \end{aligned} \right. \end{aligned}$$

(36)

where the set $\{Z_{\text {max}}, V_{\text {max}}, U_{\text {max}}\}\subset \mathbbm {R}^+$ denotes the specific values for each defined limit. The full OCP for wave energy systems can be then posed, in terms of both (35) and (36), as follows: Find an optimal control law $f^{\text {opt}}_{\text {u}}:\varXi \rightarrow \mathbbm {R}$ such that

$$\begin{aligned} \begin{aligned}&f^{\text {opt}}_{\text {u}} = \arg \max _{u\in \mathscr {U}} {\mathcal {J}}(u),\\&\text {subject to:}\\&\text {WEC nonlinear dynamics }{\tilde{\varSigma }}_{W} \text { in (31)},\\&\text {State and input constraint set }\mathscr {C} \text { in } {(36)}. \end{aligned} \end{aligned}$$

(37)

The OCP (37) is, naturally, posed in an infinite dimensional function space, which motivates the consideration of suitable tools to find an approximate solution. We consider, in this paper, a moment-based approach to compute an approximate solution of (37), by exploiting the system-theoretic concept of a moment (see e.g. [2, 3, 28]). Moments are mathematical objects which, under certain assumptions, provide a convenient parameterisation of the steady-state behaviour of system (37), for a given class of inputs, including those characterising the WEC energy harvesting process. Such a parameterisation can be explicitly used to transcribe the OCP (37) into a finite-dimensional nonlinear program (NP), i.e. in a direct optimal control fashion (see e.g. [65, 68]), as detailed in the following sections.

4.1 Direct transcription via moments

To be precise, within a moment-based approach, both external inputs, $\eta$ and $f_{\text {u}}$, are described in implicit form (see e.g [2, 27, 28]). In particular, we describe the external (uncontrollable) input affecting (31), i.e. the free-surface elevation $\eta$, in terms of a finite dimensional signal generator⁵:

(38)

with initial condition $\xi (0)$, and where $\xi (t)\in \mathbbm {R}^{\nu }$, $S\in \mathbbm {R}^{\nu \times \nu }$, and $L_{\eta }^{\intercal }\in \mathbbm {R}^{\nu }$, with $\nu = 2d$, $d>0$ an integer. The constant $\omega _0 = 2\pi /T$, with T as in (35), denotes the so-called fundamental frequency of (38) (and, correspondingly, of the free-surface elevation $\eta$). Note that Eq. (38) unveils one of the fundamental differences between the approach presented in this paper, and that in [27]; given its role in the definition of the nonlinear map g in (31), the implicit description of (38) directly takes into account the free-surface elevation $\eta$, while [27], which assumes fully linear FK effects, utilises a signal generator to describe the so-called excitation force (see Remark 14), hence ignoring any potential nonlinear coupling between $\eta$ and x in (31).

Remark 16

We assume, throughout this paper, that complete knowledge of $\eta$ is available for $t\in \varXi$, i.e. the output vector $L_{\eta }$ in (38) is known, for any given initial condition $\xi (0)$. Note that this is done without any loss of generality, and simply with the objective of avoiding the introduction of an excessively complex notation throughout the reminder of this paper. As a matter of fact, an estimate of $\eta$ can be straightforwardly incorporated into the presented framework, by following an analogous procedure to that derived in the receding-horizon moment-based control framework presented in [23].

Remark 17

The implicit description of (38) is effectively consistent with standard theory in water waves, i.e. the image of the free-surface elevation, $\eta (t)$, can be described consistently (in a statistical sense, see [52]) as the sum of $\nu$ harmonics of a fundamental frequency $\omega _0$, for any $\nu$ and T sufficiently large. If T does not fulfil this last property, such as in the case of the practical receding-horizon implementation of moment-based control (see [23]), apodisation (i.e. windowing [64]) techniques can be applied to $\eta$ within the set $\varXi$ to guarantee the well-posedness of a T-periodic extension of $\eta$, and hence the implicit form of (38) can be adopted without any loss of generality⁶. This is further discussed, from an implementation perspective, in Sect. 5.

Remark 18

The selection of a real block-diagonal form description of (38) is done without any loss of generality, and merely for representational convenience. Note that the same class of functions, i.e. state trajectories $\{\xi _i(t)\}^{\nu }_{i=1}$, can be generated with any matrix S such that $\lambda (S) = \{\pm p\omega _0\}_{p=1}^{d}\subset \mathbbm {C}^0$.

As has been demonstrated and discussed in [27], and even though the external uncontrollable input $\eta$ can be effectively expressed as in (38), the control force (to be designed via (37)) can be potentially composed of a higher number of harmonics $d + {\tilde{d}}$ in the general solution case. As a matter of fact, this ‘augmented’ representation for the control input plays a fundamental role in the quality of the approximation procedure adopted in this paper, as further discussed throughout the reminder of this section. Such a situation can be formalised via the so-called extended signal generator:

(39)

with $\overline{S}\in \mathbbm {R}^{\iota \times \iota }$, $\iota = 2(d+{\tilde{d}})$, and where the set of vectors $\{\overline{L_{\text {u}}}^{\intercal }, \overline{L_{\eta }}^{\intercal }, \overline{\xi }(t)\}\subset \mathbbm {R}^{\iota }$. Note that, given the nature of the matrix $\overline{S}$ in (39), the output vector $\overline{L}_{\eta }$ is merely the result of an appropriate inclusion, i.e.

$$\begin{aligned} \eta = L_{\eta }\xi = \left[ L_{\eta }\,0\right] \overline{\xi } = \overline{L_{\eta }}\overline{\xi }, \end{aligned}$$

(40)

with the specific inclusion map being ${\mathcal {I}}:\mathbbm {R}^{1\times \nu }\hookrightarrow \mathbbm {R}^{1\times \iota }$, ${\mathcal {I}}(L_{\eta }) \mapsto \overline{L_{\eta }}$. Furthermore, the initial condition of (39) can be defined in terms of $\xi (0)$ simply as $\overline{\xi }(0)^{\intercal } = [\xi (0)^{\intercal }\, \xi _{\text {e}}^{\intercal }(0)]$, with any vector $\xi _{\text {e}}(0)\in \mathbbm {R}^{2{\tilde{d}}}$.

Suppose, from now on (and without any loss of generality), that the initial condition of the extended signal generator (39) is set to $\overline{\xi }(0) = \varepsilon _{\iota }$, so that the pair of matrices $(\overline{S}, \overline{\xi }(0))$ is effectively excitable⁷. Then, given the nature of the WEC system (31), there exists a locally well-defined mapping $\pi :\mathbbm {R}^{\iota }\rightarrow \mathbbm {R}^2$, which solves the nonlinear partial differential equation

$$\begin{aligned} \frac{\partial \pi }{\partial \overline{\xi }} \overline{S}\overline{\xi } = f(\pi (\overline{\xi }), \overline{L_{\eta }}\overline{\xi }, \overline{L_{\text {u}}}\overline{\xi }), \end{aligned}$$

(41)

and the steady-state response map of ${\tilde{\varSigma }}_{\text {W}}$, driven by the signal generator of (39), is $x_{\text {ss}} = \pi \circ \overline{\xi }$, for any fixed state-trajectory $\overline{\xi }(t)$. We now have the tools to provide a formal definition of a moment, following, for example, [2]; from now on, we refer to the mapping $h \circ \pi$ as the moment of system (31) at the signal generator (39).

Remark 19

Local existence and uniqueness of $\pi$ in (41) can be proved via centre manifold theory (see [44, Chapter 8]), and is guaranteed both via the Poisson stable nature of the extended signal generator (39), and the Lyapunov stability properties of ${\tilde{\varSigma }}_{\text {W}}$ in (31) (see Remark 15). The interested reader is referred to, for instance, [44, Chapter 8] and [45], for a thorough discussion on the well-posedness conditions for the partial differential equation of (41).

Remark 20

For any fixed trajectory $\overline{\xi }(t)$, the moment of system (31) at the signal generator (39) coincides with the steady-state output response of ${\tilde{\varSigma }}_{\text {W}}$ driven by (39), i.e. $y_{\text {ss}}(t) = (h\circ \pi )(\overline{\xi }(t))$.

It should hopefully be clear at this point that, within the assumptions and definitions adopted within this paper, the definition of moments, and the steady-state output response of system (31), are in a one-to-one relation. In particular, the mapping $h\circ \pi$ provides a complete account of the steady-state output of ${\tilde{\varSigma }}_{\text {W}}$, for any given $\eta (t)$ and $f_{\text {u}}(t)$. We exploit such a parameterisation to approximate the solution of the OCP (37) in terms of the steady-state output map of the WEC system (31), as detailed in the following paragraphs.

We begin by noting that, though (local) existence and uniqueness of the mapping $\pi$, the solution of (41), is indeed guaranteed under the stated framework (see Remark 19), the computation of an analytical (closed-form) solution is far from trivial, if at all possible. In -other words, to make practical use of the parameterisation of the steady-state response of (31), via Remark 20, a tailored approximation method is required to compute the corresponding moment $h\circ \pi$. In the light of this, we consider the approximation procedure proposed in [27], and extend such a technique to effectively include the (input-dependent) nonlinear FK effects in the computation of the corresponding approximating moment.

In particular, [27] defines an approximation of $h\circ \pi$ by directly exploiting the nature of the extended signal generator of (39); for any fixed trajectory $\overline{\xi }(t)$, and set of output vectors $\{\overline{L_{\eta }}, \overline{L_{\text {u}}}\}$, the moment of system (31) can be approximated in terms of the solution of (39) as

(42)

where $\overline{Y} = C_{{\dot{z}}}\overline{\varPi }$ (since $h \equiv C_{{\dot{z}}}$ in the output map of (31)), with $\overline{\varPi }\in \mathbbm {R}^{2\times \iota }$ the solution of the set of algebraic equations

(43)

with ${\mathcal {T}}_{\iota } = \{t_{i}\}_{i=1}^{\iota } \subset \varXi$ a set of uniformly-distributed collocation points, and where the residual map ${\mathcal {R}}$ in (43) is defined as

(44)

Remark 21

Unlike the case posed in [27], where both the control force and so-called excitation force are regarded as a linearly combined input to the WEC system, the computation of the approximating moment $\overline{Y}$ via (43) now depends upon the free-surface elevation $\eta$ in an explicit form, i.e. two separate external inputs are being considered to solve for $\overline{Y}$.

Even though the map f, characterising the dynamics of the WEC system under nonlinear FK effects in (31), is of an integro-differential class, the system of equations (43) can be conveniently re-expressed in matrix form by virtue of moment-based theory, as discussed in the following. In particular, we begin by noting that, since $\{k_{\text {r}},k_{\text {d}}, k_{\text {FK}}^{\text {dyn}}\}\subset L^2(\mathbbm {R})$, the convolution operations in (31) can be written [21, 25], in steady-state, as

(45)

where the set of (constant) matrices $\{\overline{\mathscr {M}_{\text {r}}}, \overline{\mathscr {M}_{\text {d}}},\overline{\mathscr {M}^{\text {dyn}}_{\text {FK}}}\} \subset \mathbbm {R}^{\iota \times \iota }$, which depends upon the spectra $\lambda (\overline{S})$, can be defined as

(46)

where the mappings $K_{\text {r}}:\mathbbm {R}\rightarrow \mathbbm {C}$, $K_{\text {d}}:\mathbbm {R}\rightarrow \mathbbm {C}$, and $K_{\text {FK}}^{\text {dyn}}:\mathbbm {R}\rightarrow \mathbbm {C}$, represent the (well-defined) Fourier transform of $k_{\text {r}}$, $k_{\text {d}}$, and $k_{\text {FK}}^{\text {dyn}}$, respectively.

Furthermore, and since, for any $t\in \varXi$, the mapping f in (31) is continuous, i.e. $f\in {\mathcal {C}}^0$, then the inner-product operation between f and $\delta$ in $\varXi$ is such that $\langle f, \delta (t-t_i)\rangle = f(t_i)$. This allows for the definition of the following set of matrices:

$$\begin{aligned} \begin{aligned}&\varOmega ^{\intercal } = \begin{bmatrix} \overline{\xi }(t_1)^{\intercal } \\ \vdots \\ \overline{\xi }(t_{\iota })^{\intercal } \end{bmatrix}^{\intercal },\\&G(\overline{\varPi }, \overline{L_{\eta }})^{\intercal } = \begin{bmatrix} g(\overline{\varPi }\overline{\xi }(t_1), \overline{L_{\eta }}\overline{\xi }(t_1))^{\intercal } \\ \vdots \\ g(\overline{\varPi }\overline{\xi }(t_\iota ),\overline{L_{\eta }}\overline{\xi }(t_\iota ))^{\intercal }\\ \end{bmatrix}^{\intercal }, \end{aligned} \end{aligned}$$

(47)

where $G(\overline{\varPi }, \overline{L_{\eta }})\in \mathbbm {R}^{2\times \iota }$ and $\varOmega \in \mathbbm {R}^{\iota \times \iota }$. Finally, with the matrices introduced in (47), and the definition of the steady-state convolution ‘equivalents’ derived in (46), the set of algebraic equations posed in (43) can be written in matrix form, i.e.

$$\begin{aligned} \begin{aligned}&\overline{\varPi }\overline{S} - A\overline{\varPi } + BC_{{\dot{z}}}\overline{\varPi }\overline{\mathscr {M}_{\text {r}}} -B\overline{L_{\eta }} \left( \overline{\mathscr {M}_{\text {d}}} + \overline{\mathscr {M}_{\text {FK}}^{\text {dyn}}}\right) \\&+ B\overline{L_{\text {u}}} - G(\overline{\varPi }, \overline{L_{\eta }})\varOmega ^{-1} = 0. \end{aligned} \end{aligned}$$

(48)

Remark 22

The existence of a set of collocation points ${\mathcal {T}}_{\iota }$ such that $0\ne \lambda (\varOmega )$, i.e. $\varOmega ^{-1}$ is well-defined, is always guaranteed via the excitability property of the pair $(\overline{S},\overline{\xi }(0))$ (see [54]).

Remark 23

If the set of collocation points (time instants) ${\mathcal {T}}_{\iota }$ is chosen such that $\{ t_k\in {\mathcal {T}}_{\iota }\,\vert \, t_k = -T/2 + Tk/\iota$, $\forall k\in {\mathbb {N}}_{\iota }\}$, then the collocation approach, adopted in this paper to derive (48), identically coincides with the so-called Galerkin method [9, 30]. This, combined with the locally exponential stability of the WEC system ${\tilde{\varSigma }}_{\text {W}}$ (see Remark 15), has the following set of implications: (a) The system of equations (48) always admits a solution for any $\iota$ sufficiently large, and (b) the approximated moment converges uniformly as the number of components describing the extended signal generator (39) increases, i.e. $\overline{(h\circ \pi )}(\overline{\xi }) = \overline{Y}\overline{\xi } \rightarrow (h\circ \pi )(\overline{\xi })$ as $\iota \rightarrow + \infty$ (see [76]).

Going one step further with the analysis of the system of equations (48), and given the intrinsic relation between the state-variables z (displacement) and ${\dot{z}}$ (velocity) in (31), the algebraic equation (48) can be written in terms of $\overline{Y}$ (see “Appendix 1” for a full derivation), i.e.

$$\begin{aligned} \overline{Y}^{\intercal } + \varGamma _{\eta }\overline{L_{\eta }}^{\intercal } - \varGamma _{\text {u}}\overline{L_{\text {u}}}^{\intercal } + \varGamma _{G}\left[ G(\overline{Y}, \overline{L_{\eta }})\right] = 0, \end{aligned}$$

(49)

where the set of matrices $\{\varGamma _{\eta },\varGamma _{\text {u}}, \varGamma _{G}\}\subset \mathbbm {R}^{\iota \times \iota }$ are defined as

$$\begin{aligned} \begin{aligned} \varGamma _{\eta }&= ({\mathbb {I}}_{\iota }\otimes C_{{\dot{z}}})\varGamma _{\text {r}}^{-1}(\overline{{\mathcal {M}}_{\text {d}}}^{\intercal } + \overline{{\mathcal {M}}_{\text {FK}}^{\text {dyn}}}^{\intercal } \otimes B), \\ \varGamma _{\text {u}}&= ({\mathbb {I}}_{\iota }\otimes C_{{\dot{z}}})\varGamma _{\text {r}}^{-1}({\mathbb {I}}_{\iota }\otimes B), \\ \varGamma _{\text {G}}&= ({\mathbb {I}}_{\iota }\otimes C_{{\dot{z}}})\varGamma _{\text {r}}^{-1}(\varOmega ^{-1^{\intercal }} \otimes {\mathbb {I}}_{2}),\\ \varGamma _{\text {r}}&= (\overline{S} \otimes {\mathbb {I}}_{2}) + ({\mathbb {I}}_{\iota } \otimes A) - (\overline{{\mathcal {M}}_{\text {r}}}^{\intercal }\otimes BC_{{\dot{z}}}), \end{aligned} \end{aligned}$$

(50)

with $\varGamma _{\text {r}}\in \mathbbm {R}^{2\iota \times 2 \iota }$.

Remark 24

The following equivalent condition,

$$\begin{aligned} \overline{L_{\text {u}}}^{\intercal }= & {} \varGamma _{\text {u}}^{-1}\overline{Y}^{\intercal } + \varGamma _{\text {u}}^{-1}\varGamma _{\eta }\overline{L_{\eta }}^{\intercal } \nonumber \\&+\varGamma _{\text {u}}^{-1} \varGamma _{G}\left[ G(\overline{Y}, \overline{L_{\eta }})\right] , \end{aligned}$$

(51)

follows trivially from Eq. (49), since $0\notin \lambda (\varGamma _{\text {u}})$, due to the (local) internal stability of system (31) (see the arguments posed in [25, 26]).

As discussed at the beginning of Sect. 4, the approximated moment $\overline{(h\circ \pi )}(\overline{\xi }) = \overline{Y}\overline{\xi }$, computed via Eq. (49) (or, equivalently, Eq. (51)), can be effectively used to approximate the solution of the energy maximising OCP (37), for wave energy systems under nonlinear FK effects, in the following sense. Given the inherent one-to-one relation between moments and the steady-state output response of (31), the OCP (37) can be mapped (momentarily without considering the set $\mathscr {C}$—see Sect. 4.2 for its effective incorporation), in steady-state, to the following moment-based finite dimensional NP: For a given free-surface elevation $\eta = \overline{L_{\eta }}\overline{\xi }$, solve

$$\begin{aligned} \begin{aligned}&\overline{L_{\text {u}}}^{\text {opt}} = \arg \max _{\overline{L_{\text {u}}}^{\intercal } \in \mathbbm {R}^{\iota }} \dfrac{1}{T} \int _{\varXi } \overline{Y}\overline{\xi }(t) \overline{L_{\text {u}}}\overline{\xi }(t) \mathrm{d}t, \\&\text {subject to:}\\&\overline{Y}^{\intercal } + \varGamma _{\eta }\overline{L_{\eta }}^{\intercal } - \varGamma _{\text {u}}\overline{L_{\text {u}}}^{\intercal } + \varGamma _{G}\left[ G(\overline{Y}, \overline{L_{\eta }})\right] = 0, \end{aligned} \end{aligned}$$

(52)

where the effective optimal control input is $f_{\text {u}}^{\text {opt}} = \overline{L_{\text {u}}}^{\text {opt}}\overline{\xi }$, with $\overline{\xi }$ solution of (39).

Remark 25

Equation (52) arises by setting $y \mapsto \overline{Y}\overline{\xi }$ and $f_{\text {u}} \mapsto \overline{L_{\text {u}}}\overline{\xi }$ in (37), and by replacing the corresponding (dynamic) integro-differential WEC equality constraint (31) with the algebraic equation (49), which describes the steady-state output behaviour of WEC system (31) in terms of the corresponding approximated moment. Furthermore, note that, if $\iota \rightarrow \infty$, then the algebraic equality constraint in (52) exactly represents the steady-state output response of (31) driven by (39) (see also Remark 23).

The (now transcribed) finite-dimensional NP (52), which is a mixed optimisation problem in $\{\overline{Y},\overline{L_{\text {u}}}\}$, can be further simplified by noting that, for any admissible pair $(\overline{Y}, \overline{L_{\text {u}}})$, the following equality condition

$$\begin{aligned} \begin{aligned} \int _{\varXi } \overline{Y}\overline{\xi } \overline{L_{\text {u}}}\overline{\xi } \mathrm{d}t&= \overline{Y}\left( \int _{\varXi } \overline{\xi } \overline{\xi }^{\intercal } \mathrm{d}t\right) \overline{L_{\text {u}}}^{\intercal } = \dfrac{T}{2}\overline{Y}\overline{L_{\text {u}}}^{\intercal }, \end{aligned} \end{aligned}$$

(53)

holds. Equation (53), together with the equality derived in Remark 24, directly implies that the optimal moment-based control force $f_{\text {u}}^{\text {opt}} = \overline{L_{\text {u}}}^{\text {opt}}\overline{\xi }$ in (52) can be equivalently computed as

$$\begin{aligned} \begin{aligned} \overline{L_{\text {u}}}^{\text {opt}^{\intercal }} = \varGamma _{\text {u}}^{-1} \Big (&\overline{Y}^{\text {opt}^{\intercal }} + \varGamma _{\eta }\overline{L_{\eta }}^{\intercal } \\&+ \varGamma _{G}\left[ G(\overline{Y}^{\text {opt}^{\intercal }}, \overline{L_{\eta }})\right] \Big ), \end{aligned} \end{aligned}$$

(54)

with $\overline{Y}^{\text {opt}}$ the solution of the NP

$$\begin{aligned} \begin{aligned} \overline{Y}^{\text {opt}} = \arg&\max _{\overline{Y}^{\intercal }\in \mathbbm {R}^{\iota }} \dfrac{1}{2}\overline{Y}\varGamma _{\text {u}}^{-1} \overline{Y}^{\intercal }+\dfrac{1}{2}\overline{Y}\varGamma _{\text {u}}^{-1}\varGamma _{\eta } \overline{L_{\eta }}^{\intercal } \\&+ \dfrac{1}{2}\varGamma _{\text {u}}^{-1}\varGamma _{G}\left[ G(\overline{Y}, \overline{L_{\eta }})\right] . \end{aligned} \end{aligned}$$

(55)

We now explicitly provide a set of important remarks, which describe and discuss the nature of the moment-based solution proposed in (54)-(55).

Remark 26

For any given $\overline{L_{\eta }}$, the optimisation procedure described in (54)-(55) is carried out over the approximated moment $\overline{Y}$ only, and ‘translated’ to the effective optimal control force $f_{\text {u}}^{\text {opt}}$ via the well-defined algebraic relation (54) (see also Remark 24)⁸.

Remark 27

As in the case discussed in [27], the resulting moment-based NP in (54) is constructed as the sum of a quadratic function, and a nonlinear ‘perturbation’ term, which explicitly depends on the mapping g in (32). Note that this is indeed a highly desired representation, since the existence of a global energy-maximising control solution can be guaranteed under mild assumptions (see Remarks 28 and 29 below). In other words, with the approximation framework for nonlinear FK effects, and the subsequent control-oriented model proposed in Sect. 3, we can retain the benefits and properties of the nonlinear moment-based strategy in [27], even under the presence of nonlinear FK forces.

Remark 28

The quadratic term in (54), characterised by the Hessian matrix $\mathscr {H} = \varGamma _{\text {u}}^{-1} + \varGamma _{\text {u}}^{-1^{\intercal }}$, coincides with that derived in [25, 27], and is always strictly concave, i.e. $\lambda (\mathscr {H}) \subset \mathbbm {R}^{-}/0$. We do note that, unlike [27], which does not consider the potential existence of nonlinear FK effects, the nonlinear term in (54) now explicitly depends upon the implicit form description of the external uncontrollable input, i.e. the free-surface elevation $\eta = \overline{L_{\eta }}\overline{\xi }$.

Remark 29

Given the concave nature of the quadratic term in (55), if, for any admissible $\overline{Y}$ and $\overline{L_{\eta }}$, the nonlinear map $(\overline{Y}, \overline{L_{\eta }}) \mapsto G(\overline{Y}, \overline{L_{\eta }})$ is bounded, then the NP defined in (55) always admits a globally optimal solution⁹. This, naturally, allows the utilisation of efficient numerical optimisation routines to compute a solution for the energy-maximising optimal control law (54).

4.2 State and input constraints

The set of state and input constraints $\mathscr {C}$, defined in Eq. (36), can be incorporated into the optimisation procedure in (55) by pursuing a collocation approach. In particular, by choosing an appropriate set ${\mathcal {T}}_{\mathscr {C}} = \{t_i\}_{i=1}^{N_{\mathscr {C}}}\subset \varXi$, where the cardinality $N_{\mathscr {C}}$ is a design parameter (see also Sect. 5), the following map (see “Appendix 2” for a full derivation)

$$\begin{aligned} \mathscr {C} \mapsto \left\{ \begin{aligned} \begin{bmatrix} {\mathcal {A}}^{\intercal }_{z} \\ {\mathcal {A}}^{\intercal }_{{\dot{z}}} \\ {\mathcal {A}}^{\intercal }_{\text {u}} \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ {\mathcal {G}}_{\text {u}}(\overline{Y})^{\intercal } \end{bmatrix} \overline{Y}^{\intercal } \le \begin{bmatrix} {\mathcal {B}}^{\intercal }_{z}\\ {\mathcal {B}}^{\intercal }_{{\dot{z}}}\\ {\mathcal {B}}^{\intercal }_{\text {u}} \end{bmatrix},\\ \end{aligned} \right. \end{aligned}$$

(56)

with

$$\begin{aligned} \begin{aligned} {\mathcal {A}}_{z}&= \overline{\Delta }^{\intercal }\overline{S}^{-1^{\intercal }},\\ {\mathcal {A}}_{{\dot{z}}}&= \overline{\Delta }^{\intercal },\\ {\mathcal {A}}_{u}&= \overline{\Delta }^{\intercal }\varGamma _{\text {u}}^{-1},\\ {\mathcal {B}}_{z}&= Z_{\text {max}} \mathbf {1}_{N_{\mathscr {C}}},\\ {\mathcal {B}}_{{\dot{z}}}&= V_{\text {max}} \mathbf {1}_{N_{\mathscr {C}}},\\ {\mathcal {B}}_{\text {u}}&= U_{\text {max}} \mathbf {1}_{N_{\mathscr {C}}} - {\mathcal {A}}_{\text {u}}\varGamma _{\eta }\overline{L_{\eta }}^{\intercal },\\ {\mathcal {G}}_{\text {u}}(\overline{Y})&= {\mathcal {A}}_{\text {u}} \varGamma _G \left[ G(\overline{Y},\overline{L_{\eta }})\right] ,\\ \overline{\varLambda }&= \begin{bmatrix} \overline{\xi }(t_1)&\cdots&\overline{\xi }(t_{N_{\mathscr {C}}}) \end{bmatrix}, \\ \overline{\Delta }&= \begin{bmatrix} \overline{\varLambda }&\,&-\overline{\varLambda } \end{bmatrix}, \end{aligned} \end{aligned}$$

(57)

can be directly incorporated into (55), to explicitly take into account the set $\mathscr {C}$ in (36) within the adopted optimal control framework.

Remark 30

The pairs of matrices $({\mathcal {A}}_{z}, {\mathcal {B}}_{z})$ and $({\mathcal {A}}_{{\dot{z}}}, {\mathcal {B}}_{{\dot{z}}})$ characterise the state constraints in (36), specifically those related to displacement and velocity, respectively. Note that, within the adopted framework, such constraints are linear in the optimisation variable $\overline{Y}$, which is highly convenient from a computational efficiency perspective (see e.g. [8, 10]). In contrast, the inequality relation characterising the input (control) constraint is composed of a linear contribution, characterised by the pair $({\mathcal {A}}_{\text {u}}, {\mathcal {B}}_{\text {u}})$, and a nonlinear map $\overline{Y} \mapsto {\mathcal {G}}_{\text {u}}(\overline{Y})$, consistent with the nature of the optimisation problem (55).

5 Case study

This section presents a case study to illustrate the results and propositions presented in Sects. 3 and 4, in an integrated fashion; we first develop a data-based control-oriented model of a WEC system subject to nonlinear FK effects, and subsequently use such a model for moment-based optimal control design, aiming to effectively maximise energy extraction. In particular, to fully illustrate the features of the proposed framework, we consider a spherical heaving point absorber WEC system (see Figure 4), as extensively studied within nonlinear FK academic research, see e.g. [36, 38]. The selection of such a geometry is not only motivated by its non-uniform cross-sectional area, which clearly emphasises the relevance of nonlinear FK effects within the modelling procedure, but also by the existence of an analytical solution for the nonlinear static FK force, i.e. the map $f^{\text {st}}_{\text {FK}}$, hence providing an exact benchmark to compare the procedure, proposed in Sect. 3.2, against. Figure 4 also provides a linear frequency-domain characterisation of both radiation (top) and diffraction (bottom) impulse response functions associated with the considered heaving point absorber device, by means of a corresponding Bode plot, computed using the open-source BEM solver Nemoh [5].

Fig. 4

Schematic illustration of the spherical heaving point absorber WEC system (right), and its corresponding radiation and diffraction frequency-domain characteristics (left)

Regarding specific sea-state conditions, we consider that the WEC system in Figure 4 is subject to irregular waves stochastically characterised via JONSWAP spectra (see [41]). In particular, we consider sea-states with a significant wave height $H_w = 2.5$ [m], typical peak period $T_{w}\in [2,4]$ [s] (which fully covers the main operational range for this specific device), and a fixed peak-enhancement parameter of $\gamma = 3.3$. From now on, whenever an irregular sea-state is considered for performance assessment, the length (in seconds) associated with the time-trace of the (randomly generated) free-surface elevation, i.e. $\eta (t)$, is set to 1600 [s]. Note that such a value corresponds with more than 400 typical wave periods for each possible $T_w\in [2,4]$ [s], hence guaranteeing statistically consistent results for each considered operating condition.

The remainder of this section is organised as follows. Firstly, Sect. 5.1 presents the application of the data-based control-oriented modelling framework, introduced in Sect. 3, for both static, and dynamic nonlinear FK effects. Subsequently, Sect. 5.2 discusses the moment-based control design procedure, presented in Sect. 4, based upon the computed control-oriented WEC system, including a detailed performance assessment of such a proposed control framework.

5.1 Data-based control-oriented modelling

Recall that the control-oriented modelling framework, proposed in Sect. 3, is based exclusively on input-output data, generated by a suitable numerical nonlinear FK solver. In this paper, we consider the open-source toolbox Nlfk4all [32, 34], which provides an accurate, yet efficient (from a numerical standpoint), solution methodology for both static, and dynamic, FK effects. Different inputs are supplied to the software, according to each of the data-based modelling procedures discussed in Sect. 3, to generate a corresponding representative set of outputs, explicitly required for nonlinear FK identification.

We begin this section by describing the approximation of static FK effects, as proposed in Sect. 3.2, and, in particular, Algorithm 1. For the corresponding approximation space, i.e. the set of functions $\mathscr {P}=\{\phi _j\}$ characterising ${\tilde{f}}^{\text {st}}_{\text {FK}}$ in (13), we choose the trial set

$$\begin{aligned} \begin{aligned} \{\phi _i(z,\eta )\}^{16}_{i=1} = \{&z, z^2, z^3,\eta , \eta z, \eta z^2, \eta z^3, \eta ^2, \eta ^2 z, \\&\eta ^2 z^2, \eta ^2 z^3, \eta ^3 z, \eta ^3 z^2, \eta ^3 z^3\}, \end{aligned} \end{aligned}$$

(58)

which corresponds with the first 16 terms of the polynomial series expansion of $f^{\text {st}}_{\text {FK}}$ about $(z,\eta ) = (0,0)$. With respect to the test free-surface elevation $\eta$, used as persistently exciting input to the nonlinear FK solver to produce the corresponding set of representative identification data, we select the multisine signal $f_{\text {id}}$ specifically illustrated in Figure 1 (see also Sect. 3.1). Note that, as discussed throughout Sect. 3.2, given the static nature of system $\varSigma ^{\text {st}}_{\text {FK}}$, a single multisine test input is sufficient to characterise the proposed approximation procedure.

With respect to the specifics of the iterative procedure in Algorithm 1, we choose the sampling period, characterising the (uniformly-spaced) set ${\mathcal {T}}_{k}^w$, as 0.1 [s] (which is well within the Nyquist-Shannon limit for the multisine input signal considered), while the values for $\{w, k^0\}$ are set as $w = k^0 = 50$, which correspond to a 5 [s] window for each iteration of the proposed algorithm. The initialisation vector ${\mathcal {P}}^0\in \mathbbm {R}^{16}$, used as a starting point for Algorithm 1, is chosen randomly, aiming to highlight the convergence capabilities of the algorithm (see also Remark 11).

Fig. 5

Evolution of Algorithm 1 for each corresponding iteration (left), and final computed solution after convergence (right). Note that the obtained solution exactly coincides with that derived analytically in [36]

Fig. 6

Time-traces associated with the test (multisine) input, and resulting displacement (right), together with the corresponding value for the static FK force, plotted alongside the computed manifold

Figure 5 (left) illustrates the evolution of Algorithm 1 at each iteration k of the procedure, showing the value of each coefficient characterising the approximated mapping ${\tilde{f}}^{\text {st}}_{\text {FK}}$, according to the function space defined via (58). Note that the algorithm effectively starts the iterative procedure with random values, and converges to a solution ${\mathcal {P}}$, explicitly shown in Figure 5 (right), after $\approx$ 70 iterations. As a matter of fact, the algorithm converges to the exact analytical solution for the nonlinear static FK force effect characterising the heaving point absorber WEC of Figure 4, as derived in [36].

To finalise the presentation of the approximation results regarding static nonlinear FK effects, and explicitly illustrate the nature of the computed solution, Figure 6 shows the (sampled) time traces associated with the test (multisine) input $\eta$ (right, top), and resulting displacement z (right, bottom), together with the corresponding value for the static FK force $f^{\text {st}}_{\text {FK}}$ for each pair of samples $(\eta ,z)$, plotted alongside the manifold computed via Algorithm 1.

Fig. 7

Set of trial inputs $\mathscr {U}$ used for the control-oriented modelling of dynamic FK effects

Having presented the approximation of static nonlinear FK forces via the proposed framework, we proceed with the corresponding control-oriented modelling of dynamic FK effects, following the procedure detailed throughout Sect. 3.3, i.e. Algorithm 2. We begin such a description by providing an explicit account of the input set $\mathscr {U}$, i.e. the set of trial input signals to numerically produce a set of representative dynamic FK outputs, using the corresponding nonlinear FK solver. Recall that, within the defined operating conditions, the device is subject to stochastic wave inputs, i.e. irregular sea-states, with $H_w = 2.5$ [m] and $T_w\in [2,4]$ [s]. Using this information to construct the set of test inputs (to be able to characterise such wave operating conditions), we select $Q = 4$ different multisine signals, explicitly depicted in Figure 7. Each signal is specifically designed to emulate, via a suitable selection of the amplitudes and phases associated with each harmonic component (see Sect. 3.1 and, specifically, Remark 4), a ‘wave height’ of 2, 2.2, 2.5 and 2.7 metres, respectively. The time length of each $f_{\text {id}}(t)$ is set to $\approx 314$ [s], which corresponds to a fundamental frequency of 0.02 [rad/s]. The frequency band selected for the generation of each multisine signal is set to [0.5, 5] [rad/s], which is sufficient to cover the range of operating conditions for the WEC under scrutiny.

Regarding the specific parameters characterising Algorithm 2, the order (dimension) is set to ${\tilde{n}} = 6$, which, as detailed in the reminder of this paragraph, gives a good compromise between computational complexity and model accuracy. The set of (uniformly-spaced) trial time-shifts ${\mathcal {T}}_c$, used to compute an estimate of the (output) time-advance characterising the dynamic FK system, is such that ${\mathcal {T}}_{c}\in [0,1]$. As illustrated in Figure 8, which shows the approximation error of Algorithm 2 as a function of different values in ${\mathcal {T}}_c$, the optimal time-advance is found to be $t_c = 0.7$ [s].

Fig. 8

Approximation error of Algorithm 2 as a function of the time-advance $t_c$

Fig. 9

Bode plot of the corresponding average ETFE (dotted grey), computed approximating model (solid black), and linear BEM model (dashed red)

The results of applying the proposed algorithm are summarised in Figure 9, which presents a Bode plot including each individual ETFE $H_i(j\omega )$, average ETFE ${\bar{H}}(j\omega )$ (dotted grey), and frequency-response mapping associated with the computed approximating model ${\tilde{\varSigma }}^{\text {dyn}}_{\text {FK}}$ (solid black), which clearly presents a good fit with respect to the target ${\bar{H}}(j\omega )$. Furthermore, aiming to highlight the difference between standard linear hydrodynamic FK representations, and the proposed approach, Figure 9 also includes the frequency-response associated with the linear dynamic FK model computed via BEM solvers (Nemoh for this particular case). Note that there is a significant difference between the linear BEM model, and the representative linear structure computed via Algorithm 2, both in terms of amplitude and phase descriptions. As a matter of fact, for the latter characteristic (i.e. phase), the BEM model does not effectively capture the time-advance, having almost a zero-phase behaviour for the frequency range characterising the device operating conditions. Such misrepresentation of the phase can potentially cause a pronounced loss in energy-maximising performance for controllers based upon dynamic FK BEM models, since accurate knowledge of the instantaneous phase of the WEC system variables is fundamental in obtaining a satisfactory control performance (see e.g. [23, 81]).

Finally, and aiming to illustrate the approximation quality when effectively combining both static, and dynamic FK approximating models, Figure 10 offers a comparison between target total FK force $y_{\text {FK}}$, computed via the nonlinear FK solver Nlkf4all, and that obtained via the control-oriented approximation framework proposed in this study, i.e. ${\tilde{y}}_{\text {FK}}$, for a particular realisation of a sea-state with $T_w = 4$ [s]. In particular, Figure 10 (top) explicitly shows the time traces corresponding with target (dotted) and approximating (solid) total FK forces, while Figure 10 (bottom) provides a measure of the approximation error, computed as $(y_{\text {FK}} - {\tilde{y}}_{\text {FK}})/\max (\vert y_{\text {FK}} \vert )$, consistently showing a satisfactory approximation performance.

Fig. 10

Time-traces (top) corresponding with target (dotted) and approximating (solid) total FK force, and associated approximation error (bottom)

Fig. 11

Schematic summary of the receding-horizon implementation procedure. The green dot denotes the current time instant

5.2 Moment-based control design and assessment

Based upon the WEC control-oriented model, computed in Sect. 5.1, we now present the corresponding optimal control design for such a system subject to nonlinear FK forces. In particular, we consider the moment-based control design procedure detailed in Sect. 4. Note that, throughout the reminder of this section, the open-source toolbox Nlfk4all is always considered as the high-fidelity simulation model, so as to provide representative performance assessment results for the corresponding controller.

Concerning the control design procedure specifics, we first note that, aiming to highlight the real-time capabilities of the proposed control solution, implementation of the controller is performed in a receding-horizon fashion (see the discussion provided in Remark 17), following the framework presented in [23]. Though we do not provide a formal discussion on the theory presented in [23], the implementation procedure is briefly summarised in the following, so as to keep this paper reasonably self-contained. In particular, Figure 11 shows a schematic diagram with an overview of the main steps underlying the receding-horizon procedure to compute the corresponding control solution. Starting with knowledge (exact or estimated) of the free-surface elevation $\eta$ for a time window of length T (where energy absorption is to be maximised—see the OCP (37)), we apply a windowing¹⁰ (apodisation) mapping so as to smoothly bring $\eta$ to zero at the boundaries, and hence the derivative of its periodic extension is sufficiently smooth. The ‘windowed’ $\eta$ now automatically admits an implicit form representation in terms of the T-periodic signal generator (38), and hence the framework presented in Sect. 4 can be directly considered, i.e. the OCP (37) can be transcribed to the finite-dimensional moment-based NP (55), subject to the set of mapped state and input constraints (56). The control solution is computed, the time window is shifted by $\Delta t_{\text {rh}}$ seconds, and the procedure is repeated accordingly.

5.2.1 Controller tuning and implementation

The time-window (horizon) length T directly defines the fundamental frequency $\omega _0 = 2\pi /T$ characterising the extended signal generator (39), which, ultimately, defines the approximation space for the optimal control solution $f^{\text {opt}}_{\text {u}}$. In particular, the larger the value of T, the smaller the value of $\omega _0$, which implies a more refined ‘frequency-step’ for the definition of the optimal control solution. The selection of T is intrinsically connected to the choice of $\iota$ in (55), i.e. the number of harmonics of $\omega _0$ considered to construct the moment-based representation of the optimal energy-maximising control law and to cover the appropriate dynamic frequency range. A large value of $\iota$ increases the quality of the control solution, although having a direct impact on the computational complexity of the associated moment-based NP, which is carried over $\mathbbm {R}^{\iota }$. In practical scenarios, both T and $\iota$ can be tuned together, in terms of a single parameter, i.e. the so-called cut-off frequency $\omega _{\text {c}} = \iota \omega _0$, defining the largest multiple of $\omega _0$ used to construct the extended signal generator (39). In particular, $\omega _{\text {c}}$ can be set to a fixed value, corresponding to the largest frequency in which the stochastic description of the set of sea-states presents significant energy components. With $\omega _{\text {c}}$ fixed, and letting $\iota = \text {ceil}(T\omega _{\text {c}}/2\pi )$, we approach the tuning procedure via exhaustive (offline) simulation, by effectively changing the time-window length T while monitoring the trade-off between the value of the optimal control objective in (55) and the associated computational demand. For the particular case study presented in this section, we set the length of the time-window to $T = 15$ [s], together with $\iota = 20$, i.e. we consider a total of 10 harmonics of the fundamental frequency $\omega _0 = 2\pi /15$ [rad/s] in the signal generator (39), and hence the corresponding NP in (55) is carried out over $\mathbbm {R}^{20}$. Note that this implies a cut-off frequency for the computation of the control solution of $\omega _{\text {c}} \approx 4$ [rad/s].

The receding time step is set to $\Delta t_{\text {rh}} = 0.1$ [s], which corresponds to an order of magnitude below the typical sampling time of a full-scale WEC device (see e.g. [23]), consistent with standard real-time requirements. The set of collocation points ${\mathcal {T}}_{\mathscr {C}}$, used to enforce the mapped constraints in (56), is tuned by selecting a uniformly distributed set of time instants, with the same time step chosen for $\Delta t_{\text {rh}}$, i.e. 0.1 [s]. For the case study presented in this section, the cardinality of ${\mathcal {T}}_{\mathscr {C}}$ is hence $N_{\mathscr {C}} = T/0.1 = 150$. We do note that, while one could consider a ‘larger’ time step for the constraint collocation instants in order to enhance computational requirements, this can potentially affect the constraint enforcement capabilities of the controller.

The values for each specific state and input limitation are defined according to the technical specifications of the device under scrutiny. For this case study, we set $Z_{\text {max}} = 2.5$ [m] (displacement constraint), which prevents the device from either being fully submerged, or fully out of the water, $V_{\text {max}} = 2.5$ [m/s] (velocity constraint), and $U_{\text {max}} = 1.5\times 10^5$ [N] (control force constraint), for all the considered sea-states.

Regarding the specific computation of the associated numerical solution, the algorithm used to solve the NP (55), subject to the set of state and input constraints in (56), is based on the interior-point method described in [80], implemented in MATLAB Simulink. As can be appreciated in Figure 12, the average computational time required to compute the control solution corresponding to each receding-horizon step (solid-green line), i.e. each time window of length T, is $\approx 10^{-2}$ [s], which is one order of magnitude smaller than $\Delta t_{\text {rh}}$ (solid-red line), hence always consistently achieving real-time performance.

Fig. 12

Computational time required by each corresponding control calculation (indicated using ‘$+$’ markers) as a function of the simulation time. The solid-green line shows the average required computational time, while the solid-red line indicates the real-time limit, defined by $\Delta t_{\text {rh}}$

5.2.2 Performance assessment

Before effectively presenting and discussing performance results, and aiming to provide a comparison of the proposed strategy against a benchmark WEC controller, we introduce a well-established control methodology in the WEC literature, i.e. the so-called ‘reactive control’, which is essentially a proportional-integral (PI) control structure $K_{\text {PI}}$ achieving impedance-matching (see e.g. [19, 74]) at a (suitably selected) single interpolation frequency. To be precise, $K_{\text {PI}}:\mathbbm {C}\rightarrow \mathbbm {C}$, $s\mapsto K_{\text {PI}}(s)$, is defined as

$$\begin{aligned} K_{\text {PI}}(s) = \theta _1 + \dfrac{\theta _2}{s}, \end{aligned}$$

(59)

where the set $\{\theta _1, \theta _2\}\subset \mathbbm {R}$ can be uniquely selected to interpolate the so-called optimal control impedance $I_{\text {u}}$, at a given interpolation frequency $\omega _{\text {I}}$ (the reader is referred to “Appendix 3” for the explicit definition of the control impedance $I_{\text {u}}$ and the corresponding interpolation procedure).

Fig. 13

Frequency-response mapping for the optimal control impedance (solid black), and each corresponding $K_{\text {PI}}$ (dotted), designed to interpolate $I_{\text {u}}$ at each considered peak period $T_{w}\in [2,4]$

Aiming to consider a consistent benchmark comparison case, we consider an optimal PI control structure (59) for each specific sea-state, i.e. we re-design the set $\{\theta _1,\theta _2\}$ as a function of the specific operating condition, instead of simply fixing a single controller $K_{\text {PI}}$ for all the considered sea-states. In particular, we design (59) to interpolate the optimal control impedance (68) at $\omega _{\text {I}} = 2\pi /T_w$, for each considered $T_w\in [2,4]$. This situation is explicitly illustrated in Figure 13, where the frequency-response of the optimal control impedance is shown (solid black), alongside the set of Bode plots characterising $K_{\text {PI}}$ (dotted), designed to interpolate $I_{\text {u}}$ for each corresponding (input) sea-state condition. Note that, in (natural) resonance conditions ($T_{w} \approx 3$ [s]), the optimal $K_{\text {PI}}$ is simply a constant value, i.e. it becomes passive.

Fig. 14

Energy absorption results for both reactive PI (red), and the proposed nonlinear moment-based controller (green)

Figure 14 displays the main performance results obtained with the proposed nonlinear optimal moment-based controller, based upon the control-oriented modelling framework presented in this study. In particular, Figure 14 shows energy absorption for both reactive PI (red), and nonlinear moment-based control (green), computed for the full time-length of each free-surface elevation (1600 [s]). Note that the proposed controller is able to significantly outperform the benchmark case (which has been optimised for each particular typical peak period), by means of a single set of control parameters, hence directly showing the capabilities of the nonlinear moment-based approach to maximise energy in a wide range of sea-states (i.e. operating conditions).

Fig. 15

Applied control force (top), displacement of the WEC device under controlled conditions (centre), and instantaneous power (bottom), for a random realisation of a sea-state with $T_w = 3.5$ [s]

Finally, and aiming to further illustrate and discuss the comparative controller performance, Figure 15 shows applied control force (top), displacement of the WEC device under controlled conditions (centre), and instantaneous power (bottom), for a random realisation of a sea-state with $T_w = 3.5$ [s]. With respect to effective control forces (top), a consistent difference in instantaneous phase can be appreciated between the optimised PI controller and the proposed nonlinear moment-based controller, which leads to increased displacement values under controlled conditions (centre), and a higher requirement (in average) of negative (reactive) instantaneous power flow (bottom). In other words, the proposed control solution is able to absorb a significantly higher value of energy, while requiring, at the same time, more conservative displacement ranges, and less reactive power flow to achieve optimality.

6 Conclusions

In this paper, an integrated framework for optimal control of wave energy systems, subject to nonlinear FK forces, is proposed, by fulfilling two fundamental objectives. We first derive a data-based control-oriented modelling scheme, capable of computing suitable models for control purposes, by means of representative input/output data, numerically obtained via an appropriate FK solver (Nlfk4all for our study). Both static and dynamic FK effects are included within such a modelling framework, by means of tailored approximation procedures. As illustrated in the case study presented in Sect. 5, the proposed algorithm for nonlinear FK forces is capable of recovering the exact analytical solution based solely on data, while the methodology outlined in this study for the dynamic FK case delivers a representative linear model which is effectively able to capture the associated dynamics accurately, unlike its BEM-based counterpart.

Secondly, we consider the computed model to design an optimal moment-based controller, capable of effectively maximising energy absorption for devices under such nonlinear hydrodynamic forces. The proposed controller is shown to outperform a benchmark strategy, well-established within the WEC control literature, in terms of energy absorption (with an increase of up to 3 times in performance), while effectively incorporating state and input constraints, and more conservative requirements in terms of operational space (i.e. motion range), while requiring less reactive (bi-directional) power flow to achieve optimality. Furthermore, by virtue of the efficient steady-state parameterisation offered by moment-based theory, the solution of the corresponding nonlinear program can be performed in real-time, hence not only providing a solid theoretical framework to achieve energy-maximisation, but also a practical control solution to support the pathway towards effective commercialisation of wave energy systems.

Appendices

Appendices

Appendix 1: On the derivation of Eq. (49)

We begin by recalling a well-known property of the vec operator (see e.g. [43]), i.e. the following equivalence

$$\begin{aligned} \left[ A_1A_2A_3\right] = (A_{3}^{\intercal }\otimes A_1) \left[ A_2\right] , \end{aligned}$$

(60)

for any $A_1\in \mathbbm {R}^{k\times l}$, $A_2\in \mathbbm {R}^{l\times m}$, and $A_3\in \mathbbm {R}^{m\times n}$, holds.

An application of the vec operator to Eq. (48), together with the property defined via (60), yields

$$\begin{aligned} \begin{aligned}&\left( -\overline{S}^{\intercal }\otimes BC_{{\dot{z}}}\right) \left[ \overline{\varPi }\right] + \left( \overline{\mathscr {M}_{\text {d}}} + \overline{\mathscr {M}_{\text {FK}}^{\text {dyn}}} \otimes B\right) \left[ \overline{L_{\eta }}\right] \\&- \left( {\mathbb {I}}_{\iota }\otimes B\right) \left[ \overline{L_{\text {u}}}\right] + \left( \varOmega ^{-1}\otimes {\mathbb {I}}_2\right) \left[ G(\overline{\varPi }, \overline{L_{\eta }})\right] . \end{aligned} \end{aligned}$$

(61)

Equation (49), together with the corresponding matrix definitions in (50), follow by multiplying (61) by $({\mathbb {I}}_{\iota }\otimes C_{{\dot{z}}})$, and noting that

$$\begin{aligned} \begin{aligned}&({\mathbb {I}}_{\iota }\otimes C_{{\dot{z}}})\left[ \overline{\varPi } \right] = \left[ C_{{\dot{z}}}\overline{\varPi }\right] = \left[ \overline{Y}\right] = \overline{Y}^{\intercal }, \\&\left[ \overline{L_{\eta }}\right] = \overline{L_{\eta }}^{\intercal },\\&\left[ \overline{L_{\text {u}}}\right] = \overline{L_{\text {u}}}^{\intercal }, \end{aligned} \end{aligned}$$

(62)

where we have used the fact that $\overline{S} = -\overline{S}^{\intercal }$. To finalise proving the claim, we note that the map $(\overline{\varPi }, \overline{L_{\eta }}) \mapsto G(\overline{\varPi }, \overline{L_{\eta }})$ can be fully characterised in terms of $\overline{Y}$, i.e. the map¹¹$(\overline{Y}, \overline{L_{\eta }}) \mapsto G(\overline{Y}, \overline{L_{\eta }})$, by virtue of the following relation. Recall that the state-vector x of system (31) is $x = [z,\, {\dot{z}}] = [z\, y]$, and that, in steady-state conditions, $x_{\text {ss}} = \overline{\varPi }\overline{\xi }$ and $y_{\text {ss}} = C_{{\dot{z}}}\overline{\varPi }\overline{\xi } = \overline{Y}\overline{\xi }$. Following [25, 69], one can derive that

$$\begin{aligned} z = \int y\,\mathrm{d}t \mapsto z_{\text {ss}} = \overline{Y}\overline{S}^{-1}\overline{\xi }, \end{aligned}$$

(63)

and hence $\overline{\varPi }$ can be fully written in terms of $\overline{Y}$ as

$$\begin{aligned} \overline{\varPi } = \begin{bmatrix} \overline{Y}\overline{S}^{-1} \\ \overline{Y} \end{bmatrix}. \end{aligned}$$

(64)

Appendix 2: On the derivation of Eq. (56)

Recall the definition of the state and input constraints $\mathscr {C}$ in (36). In steady-state conditions, one can explicitly use the approximated moment $\overline{Y}$ to map the set $\mathscr {C}$ as

$$\begin{aligned} \mathscr {C}:\left\{ \begin{aligned}&\vert z(t) \vert \le Z_{\text {max}},\\&\vert {\dot{z}}(t) \vert \le V_{\text {max}},\\&\vert u(t) \vert \le U_{\text {max}}, \end{aligned} \right. \mapsto \mathscr {C}_{\text {ss}} \left\{ \begin{aligned}&\vert \overline{Y}\overline{S}^{-1}\overline{\xi }(t) \vert \le Z_{\text {max}},\\&\vert \overline{Y}\overline{\xi }(t) \vert \le V_{\text {max}},\\&\vert \overline{L_{\text {u}}}\overline{\xi }(t) \vert \le U_{\text {max}}, \end{aligned}\right. \end{aligned}$$

(65)

where we have explicitly used the relation posed in equation (63). Using the property that, for any $f\in {\mathcal {C}}^1$ and $\alpha \in \mathbbm {R}^+/0$, $\vert f(t)\vert \le \alpha \rightarrow f(t) \in [-\alpha , \alpha ]$, and enforcing the set $\mathscr {C}_{\text {ss}}$ at the collocation set ${\mathcal {T}}_{\mathscr {C}}\subset \mathbbm {R}^+$, we can write (65) as

$$\begin{aligned} \begin{aligned}&\overline{\Delta }^{\intercal } \overline{S}^{-1^{\intercal }}\overline{Y}^{\intercal } \le Z_{\text {max}}\mathbf {1}_{N_{\mathscr {C}}},\\&\overline{\Delta }^{\intercal } \overline{Y}^{\intercal } \le V_{\text {max}}\mathbf {1}_{N_{\mathscr {C}}},\\&\overline{\Delta }^{\intercal }\overline{L_{\text {u}}}^{\intercal } \le U_{\text {max}} \mathbf {1}_{N_{\mathscr {C}}}, \end{aligned} \end{aligned}$$

(66)

from which the pairs of matrices $({\mathcal {A}}_{z}, {\mathcal {B}}_{z})$ and $({\mathcal {A}}_{{\dot{z}}}, {\mathcal {B}}_{{\dot{z}}})$ in (57) follow directly. Finally, the expressions for the pair $({\mathcal {A}}_{\text {u}}, {\mathcal {B}}_{\text {u}})$, and the nonlinear map $\overline{Y} \mapsto {\mathcal {G}}_{\text {u}}(\overline{Y})$ in (57), can be obtained by noting that

$$\begin{aligned} \overline{\Delta }^{\intercal }\overline{L_{\text {u}}}^{\intercal }= & {} \overline{\Delta }^{\intercal }\varGamma _{\text {u}}^{-1}\left( \overline{Y}^{\intercal } + \varGamma _{\eta }\overline{L_{\eta }}^{\intercal }\right. \nonumber \\&\left. + \varGamma _{G}\left[ G(\overline{Y}, \overline{L_{\eta }})\right] \right) , \end{aligned}$$

(67)

for any fixed $\overline{L_{\eta }}$, via Remark 24.

Appendix 2: On the design of (59)

The PI control structure in (59) is designed so as to interpolate the so-called optimal control impedance $I_{\text {u}}:\mathbbm {C}^0 \rightarrow \mathbbm {C}$, $j\omega \mapsto I_{\text {u}}(j\omega )$, which can be simply defined [19, 31] as

$$\begin{aligned} I_{\text {u}}(j\omega ) = \dfrac{1}{W^{\star }(j\omega )}, \end{aligned}$$

(68)

where ${}^{\star }:\mathbbm {C}\rightarrow \mathbbm {C}$ denotes the Hermitian operator, and the mapping $W:\mathbbm {C}^0\rightarrow \mathbbm {C}$ is the frequency-response of the linearised WEC model (34), i.e.

$$\begin{aligned} W(j\omega ) = C_{{\dot{z}}}\left( j\omega - A + BC_{{\dot{z}}}K_{\text {r}}(j\omega ) \right) ^{-1}B. \end{aligned}$$

(69)

Given any interpolation frequency $\omega _{\text {I}}$, the design of the set of parameters $\{\theta _1, \theta _2\}$ in (59) is performed such that the following equality

$$\begin{aligned} K_{\text {PI}}(j\omega _{\text {I}}) = I_{\text {u}}(j\omega _{\text {I}}), \end{aligned}$$

(70)

holds, which, for the specific structure in (59), can be straightforwardly achieved with the following choice of parameters:

$$\begin{aligned} \theta _1 = \Re \left( I_{\text {u}}(j\omega )\right) , \quad \theta _2 = -\omega _{\text {I}}\Im \left( I_{\text {u}}(j\omega )\right) . \end{aligned}$$

(71)