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

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.


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 closedform expression for the so-called latching time (which is ill-posed in the case of irregular wave inputs 1 ).
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 incompati-ble 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 controloriented 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-theart 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 momentbased 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 controloriented 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 welldefined 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.

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.

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

Scalars, vectors and matrices
The symbol 0 stands for any zero element, dimensioned according to the context. The symbol I n denotes the identity matrix in C n×n , while the notation 1 n×m is used to denote a n × m Hadamard identity matrix (i.e. a n × m matrix with all its entries equal to 1). The superscript denotes the transposition operator. The spectrum of a matrix A ∈ R n×n , i.e. the set of its eigenvalues, is denoted by λ(A). The symbol denotes the direct sum of n matrices, i.e. n i=1 A i = diag(A 1 , A 2 , . . . , A n ). The notation (z) and (z), with z ∈ C, stands for the real-part and the imaginarypart operators, respectively. The Kronecker product between two matrices M 1 ∈ R n×m and M 2 ∈ R p×q is denoted by M 1 ⊗ M 2 ∈ R np×mq . The vectorisation operator acting on a matrix A ∈ C n×m is denoted as vec(A) = [A] ∈ C nm . The symbol ε n ∈ R n denotes a vector with odd entries equal to 1 and even entries equal to 0.

Functions
The generalised Dirac-δ function is denoted as δ : R → R, t → δ(t). Given two functions f 1 and f 2 , such that f 1 : X → Y and f 2 : Z → X , the composition f 1 ( f 2 (z)), which maps all z ∈ Z to f 1 ( f 2 (z)) ∈ Y , is denoted as f 1 • f 2 . The convolution between two functions g 1 and g 2 over R, i.e. R g 1 (τ )g 2 (t − τ )dτ is denoted as g 1 * g 2 . Let w 1 and w 2 be two functions belonging to the space L 2 (Ξ ), with Ξ ⊂ R closed. Then, the inner-product between w 1 and w 2 is given by w 1 , w 2 = Ξ w 1 (τ )w 2 (τ )dτ .

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 degreeof-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 : R + → R, t → z(t) and η : R + → R, t → η(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 Σ W written, for t ∈ R + , [15,29,37] as where m ∈ R + is the generalised mass of the device, f r : R → R is the radiation force, f v : R → R represents viscous effects, f d : R → R is the diffraction force, and the mappings f st FK : R × R → R and f dyn FK : R × R → 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 : R + → R is assumed to be a linear combination of device displacement and velocity, defined via the matrix C ∈ R 2 .
The radiation force f r is modelled based on linear potential theory and, using the well-known Cummins' equation [13], can be written as where m ∞ is the so-called added-mass at infinite frequency, and k r : R + → R, k r ∈ L 2 (R), is the (causal) radiation impulse response function describing the memory effect of the fluid response.

Remark 2
The impulse response function k r fully characterises a linear time-invariant system with inputż and output f r . The rationality behind representing f r in terms of a convolution operator stems from the fact that the impulse response k 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 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 v , which represents viscous effects, is written in terms of a smooth approximation of the so-called Morison equation [53], i.e.
with ∈ R + sufficiently small, and α v ∈ R + directly depending on the geometric properties of the device. The diffraction force, represented via f d , can be described (analogously to the radiation force equation (2)), in terms of a convolution operator, i.e.
where the impulse response k d : R + → R, k d ∈ L 2 (R), fully characterises a linear time-invariant system with input η and output f d . Note that, as in the case of the radiation force (see Remark 2), the diffraction kernel k d is virtually always computed numerically via BEM solvers, i.e. in a non-parametric form, hence the rational behind representing f d in terms of an appropriate convolution operator.
Finally, the map f u : R + → R, t → f 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 momentbased approach, in Sect. 4.

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 socalled 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 S w (η, z) ≡ S w be the instantaneous wetted surface of the device, and let the mappings p st : R → R, and p dyn : R × R → R denote the static and dynamic pressures, respectively, defined as where Ψ I denotes the so-called incident potential function (see [48,55]). The integration of the static pressure p st , as defined in Eq. (5), over the instantaneous wetted surface S w , yields the so-called static FK force, i.e. the map f st FK , defined as where f g denotes the gravity force, and n z is a (normalised) vector orthogonal to S w , while the corresponding integration of the dynamic pressure mapping p dyn in (5), i.e.
gives rise to the so-called dynamic Froude-Krylov force f dyn FK . In the case of a generic geometry, i.e. a general surface S w , the computation of both (6) and (7) can be performed by an appropriate spatial discretisation (e.g. mesh-based) of 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 2 .
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 controloriented, 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].

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.

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 3 , due to the specific banded nature of the free-surface elevation (wave) input η. We give a formal definition of such signals in the following paragraph. A multisine f id is a periodic signal with a bandlimited spectrum. As such, it can be fully represented in terms of a Fourier series, i.e. a trigonometric sum of order K , with K ∈ N finite, such as where where ω id = 2π/T id , with T id the measurement period, and l p a positive integer. The frequency ω 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 {ψ p } is optimised aiming at minimising the socalled 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 ∈ R + , ∀ p ∈ N K ), [73] suggests the set of phases in (8) as

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 {φ p }, which are often based upon tailored optimisation routines. The reader is referred to, for instance, [40], for further discussion on this topic.
Aiming to briefly illustrate the nature of such a signal, Figure 1 (top) presents a multisine signal with measurement period T id ≈ 314 [s], which corresponds to a fundamental frequency of ω id ≈ 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 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.

Static FK effects
The mapping f st FK is a static function which depends upon both the free-surface elevation η, and the displacement of the device z. Although, for some specific geometries, f st 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 st FK can be seen, from a system-theoretic perspective, as a static system Σ st , with both η and z as inputs, i.e.
where y st (t) ∈ 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 which is essentially system Σ W in (1) without considering static FK effects, and 'decoupling' Σ 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 u = 0, ∀t ∈ R + . Note that f u does not play a role in the identification of (10), since Σ st is fully characterised by a static mapping depending only directly on η and z.
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 st FK belongs to the space generated by a family of real-valued functions for every {η(t), z(t)} ⊂ 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 st FK , as detailed in the following: From now on, we call the mappingf st FK , defined as with N ∈ 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 st 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 P = {φ j } N j=1 . 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Σ st , analogously to Eq. (10), simply as where, clearly,ỹ st ≈ y st in the sense of (13).

Remark 8
In practice, the set of functions P can be selected via a trial and error procedure, starting with, for instance, a polynomial expansion, if the mapping f st FK is known to be smooth.
In order to proceed with the proposition of the databased approximation framework for static FK forces, we introduce the following set of auxiliary variables: where {P , Φ(η, z)} ⊂ R N , so that the approximated mapping defined in (13) can be written in a compact form as With the compact definition presented in (16), the approximation problem reduces to finding a suitable matrix P, for any given basis-function vector Φ(η, z).
In particular, we now propose a method to compute P based upon a recursive least-squares approach, inspired by the data-driven model reduction framework of [70].
. . , t k−1 , t k } ⊂ R + be a set of time-instants, where we numerically evaluate the output of the target static FK system, i.e. y st in (10). Note that the set T w k basically represents a moving window of w ∈ N samples. Suppose that, in the numerical evaluation of y st , we select the external input η (i.e. freesurface elevation) as a band-limited multisine signal, i.e. η = f id (see Sect. 3.1), and let Ξ k ∈ R w×N and Υ k ∈ R w , with w ≥ N , be defined as Let P k be an online estimate of the matrix 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 P k based upon the following optimal criterion: i.e. to find the least-squares solution of Ξ k P k = Υ 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 Ξ k is full column rank for all k ∈ N.
Remark 9 Note that, since η 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 leastsquares solution to solve for 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 id for the generation of η, and a sufficiently large initial value k 0 , the matrices Ξ k and Υ k are constructed iteratively, using information on the (supplied) input η = f id , and the numerically computed device motion z, and total static force y 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 .
Algorithm 1 (Static FK approximation) Let k 0 be a sufficiently large integer, and let be a sufficiently small user-defined error tolerance. Define P 0 ∈ R 1×N as the initialisation vector for the computation of P.
Construct the matrices Ξ k and Υ k as in (17);

Remark 10
The matrix Ψ k is always well-defined due to the fact that Ξ k is full column rank for all k ∈ N.

Remark 11
The selection of 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 P are updated at each k), after a sufficient number of iterations the effect of such a selection is effectively forgotten.

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 dyn FK can be generally represented as the output mapping y dyn of a dynamical system Σ dyn , with both η and z as inputs, i.e.

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Σ 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Σ 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 blackbox structures, using only input-output data in the frequency-domain. Such a methodology is discussed in the following paragraphs. Let be a set of suitably selected multisine input signals (i.e. free-surface elevation profiles), with Q ∈ N ≥1 , generating a corresponding set denote an input-output pair of signals for system (21). We define the so-called empirical transfer function estimate (ETFE) H i : C 0 → C, jω → H i ( jω), for each input-output pair, in terms of the expression with i ∈ N Q , and where N i : C 0 → C and Y dyn i , one can readily obtain the so-called average ETFE,H , computed with the aim of building a low-variance set to use as input to the frequencydomain identification algorithm (see e.g. [63]). The explicit computation ofH is simply given by: Recall that the ultimate objective of the proposed system identification procedure is to obtain a parametric form which approximates the behaviour of Σ dyn in terms of a representative linear structureΣ 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 where the term e jωt c denotes the frequency-domain equivalent of the time-shift (advance) corresponding with the dynamic FK system (see Remark 12) and, hence,H c only represents the causal component ofH . The strategy considered here, to compute a statespace 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 inputoutput information encoded in the mappingH ) 4 . 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.
be a set of trial time-shifts, with a l sufficiently small, and a h and P ∈ N sufficiently large. Let Σ denote an approximated state-space system, computed from an average ETFEH , with a (finite) userselected order (dimension)ñ, and letH denote its associated transfer function.
with H i as in (22); Compute the average ETFEH as in (23); 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 T c is chosen. For each trial value contained in the set,H c is constructed (as in equation (24)), and a corresponding approximating state-space structureΣ of order (dimension)ñ is computed using subspace techniques. Note that the latter step is indicated in Algorithm 2 with the 'function' identify(· , ·). Finally, the systemΣ i producing the lowest fitting error E i , together with its associated t c value, is selected as the linear representative approximating modelΣ 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 FK = y st + y dyn , viaỹ FK =ỹ st +ỹ dyn .

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 controloriented modelΣ W , approximating the nonlinear WEC dynamics Σ 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 mappingf st FK , characterising nonlinear static FK forces within the presented approximation framework, can be 'separated' as wherẽ i.e. into linearf st FK l , and nonlinearf st FK nl , contributions, respectively.

Remark 13
The linear mapf st FK l depends only on z: From the physical laws governing the WEC dynamics, it is possible to show that ∂Φ ∂η (0,0) = 0, and hence the linear part of the static FK force depends only upon the device displacement. In fact, the coefficient P ∂Φ ∂z (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 C 1 function g W : R × R × R → R, defined as 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Σ dyn in terms of its associated impulse response function, i.e.
wherek dyn FK : R + → R,k dyn FK ∈ L 2 (R), can be directly defined in terms of the triple of matrices (F, G, H ) composing (21) as (29) Using the results presented up until this point, the approximating WEC dynamical equation can be conveniently written as with M = m + m ∞ . Note that the linear and nonlinear contributions are explicitly separated in equation (30). Finally, and to provide a model compatible with control design (as approached in Sect. 4), let x = [z,ż] ∈ R 2 be the state vector associated with the WEC system. The approximating control-oriented model can be then expressed as with the map g : R × R 2 → R 2 given by and where the pair of matrices (A, B), defining (31), can be straightforwardly defined as 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 [On the stability ofΣ W ] It is straightforward to show that (x, η, f u ) = (0, 0, 0) is an equilibrium point for (31). Furthermore, given that ∂g(x,η) ∂ x | (0,0) = 0, with g as in (32), the local stability properties of the zero-equilibrium ofΣ W are governed by the linearised equatioṅ 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 r (see e.g. [29,75]).

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 wavedependent 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 opt 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 J : U u → R, u → J (u), defined (see e.g. [18,67]) as where y is as defined in (31), Ξ = [0, T ] ⊂ R + , and U 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 C can be formally written as C : where the set {Z max , V max , U max } ⊂ 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) 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.

Direct transcription via moments
To be precise, within a moment-based approach, both external inputs, η and f 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 η, in terms of a finite dimensional signal generator 5 : with initial condition ξ(0), and where ξ(t) ∈ R ν , S ∈ R ν×ν , and L η ∈ R ν , with ν = 2d, d > 0 an integer. The constant ω 0 = 2π/T , with T as in (35), denotes the so-called fundamental frequency of (38) (and, correspondingly, of the free-surface elevation η). 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 η, while [27], which assumes fully linear FK effects, utilises a signal generator to describe the socalled excitation force (see Remark 14), hence ignoring any potential nonlinear coupling between η and x in (31).

Remark 16
We assume, throughout this paper, that complete knowledge of η is available for t ∈ Ξ , i.e. the output vector L η in (38) is known, for any given initial condition ξ(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 η can be straightforwardly incorporated into the presented framework, by following an analogous procedure to that derived in the recedinghorizon 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, η(t), can be described consistently (in a statistical sense, see [52]) as the sum of ν harmonics of a fundamental frequency ω 0 , for any ν 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 η within the set Ξ to guarantee the well-posedness of a T -periodic extension of η, and hence the implicit form of (38) can be adopted without any loss of generality 6 . 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 {ξ i (t)} ν i=1 , can be generated with any matrix S such that As has been demonstrated and discussed in [27], and even though the external uncontrollable input η 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 +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: with S ∈ R ι×ι , ι = 2(d +d), and where the set of vectors {L u , L η , ξ(t)} ⊂ R ι . Note that, given the nature of the matrix S in (39), the output vector L η is merely the result of an appropriate inclusion, i.e.
with the specific inclusion map being I : R 1×ν → R 1×ι , I(L η ) → L η . Furthermore, the initial condition of (39) can be defined in terms of ξ(0) simply as ξ(0) = [ξ(0) ξ e (0)], with any vector ξ e (0) ∈ R 2d . Suppose, from now on (and without any loss of generality), that the initial condition of the extended signal generator (39) is set to ξ(0) = ε ι , so that the pair of 6 The reader is referred to [23] for the receding-horizon implementation of moment-based control. matrices (S, ξ(0)) is effectively excitable 7 . Then, given the nature of the WEC system (31), there exists a locally well-defined mapping π : R ι → R 2 , which solves the nonlinear partial differential equation ∂π and the steady-state response map ofΣ W , driven by the signal generator of (39), is x ss = π • ξ , for any fixed state-trajectory ξ(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 • π as the moment of system (31) at the signal generator (39).

Remark 19
Local existence and uniqueness of π 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Σ 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).
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 • π provides a complete account of the steady-state output ofΣ W , for any given η(t) and f 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 π , the solution of (41), is indeed guaranteed under the stated framework (see Remark 19), the computation of an analytical (closedform) 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 • π . 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 • π by directly exploiting the nature of the extended signal generator of (39); for any fixed trajectory ξ(t), and set of output vectors {L η , L u }, the moment of system (31) can be approximated in terms of the solution of (39) as where Y = CżΠ (since h ≡ Cż in the output map of (31)), with Π ∈ R 2×ι the solution of the set of algebraic equations ⊂ Ξ a set of uniformly-distributed collocation points, and where the residual map R in (43) is defined as 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 Y via (43) now depends upon the free-surface elevation η in an explicit form, i.e. two separate external inputs are being considered to solve for 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 r , k d , k dyn FK } ⊂ L 2 (R), the convolution operations in (31) can be written [21,25], in steady-state, as where the set of (constant) matrices {M r , M d , M dyn FK } ⊂ R ι×ι , which depends upon the spectra λ(S), can be defined as , (46) where the mappings K r : R → C, K d : R → C, and K dyn FK : R → C, represent the (well-defined) Fourier transform of k r , k d , and k dyn FK , respectively. Furthermore, and since, for any t ∈ Ξ , the mapping f in (31) is continuous, i.e. f ∈ C 0 , then the innerproduct operation between f and δ in Ξ is such that . This allows for the definition of the following set of matrices: . . .
where G(Π, L η ) ∈ R 2×ι and Ω ∈ R ι×ι . 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.

Remark 23
If the set of collocation points (time instants) T ι is chosen such that {t k ∈ T ι | t k = −T /2 + T k/ι, ∀k ∈ N ι }, 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 sys-temΣ W (see Remark 15), has the following set of implications: (a) The system of equations (48) always admits a solution for any ι sufficiently large, and (b) the approximated moment converges uniformly as the number of components describing the extended signal generator (39) [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) anḋ z (velocity) in (31), the algebraic equation (48) can be written in terms of Y (see "Appendix 1" for a full derivation), i.e. (49) where the set of matrices {Γ η , Γ u , Γ G } ⊂ R ι×ι are defined as with Γ r ∈ R 2ι×2ι .
As discussed at the beginning of Sect. 4, the approximated moment (h • π)(ξ) = Y ξ , 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 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 η = L η ξ , solve subject to: where the effective optimal control input is f opt u = L u opt ξ , with ξ solution of (39). (52) arises by setting y → Y ξ and f u → L u ξ 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 ι → ∞, then the algebraic equality constraint in (52) exactly represents the steady-state output response of (31) driven by (39) (see also Remark 23).

Remark 25 Equation
The (now transcribed) finite-dimensional NP (52), which is a mixed optimisation problem in {Y , L u }, can be further simplified by noting that, for any admissible pair (Y , L u ), the following equality condition holds. Equation (53), together with the equality derived in Remark 24, directly implies that the optimal momentbased control force f opt u = L u opt ξ in (52) can be equivalently computed as with Y opt the solution of the NP We now explicitly provide a set of important remarks, which describe and discuss the nature of the momentbased solution proposed in (54)- (55).

Remark 26
For any given L η , the optimisation procedure described in (54)-(55) is carried out over the approximated moment Y only, and 'translated' to the effective optimal control force f opt u via the well-defined algebraic relation (54) (see also Remark 24) 8 .

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 energymaximising 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 controloriented 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 H = Γ −1 u + Γ −1 u , coincides with that derived in [25,27], and is always strictly concave, i.e. λ(H ) ⊂ 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 η = L η ξ .

Remark 29
Given the concave nature of the quadratic term in (55), if, for any admissible Y and L η , the nonlinear map (Y , L η ) → G(Y , L η ) is bounded, then the NP defined in (55) always admits a globally optimal solution 9 . This, naturally, allows the utilisation of efficient numerical optimisation routines to compute a solution for the energy-maximising optimal control law (54).

State and input constraints
The set of state and input constraints 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 where the cardinality N C is a design parameter (see also Sect. 5), the following map (see "Appendix 2" for a full derivation) 9 This statement follows from the so-called Γ -convex nature of the objective function in (54) when the map G is bounded. The interested reader is referred to, for example, [60][61][62] for further detail. with can be directly incorporated into (55), to explicitly take into account the set C in (36) within the adopted optimal control framework.

Remark 30
The pairs of matrices (A z , B z ) and (Aż, Bż) 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 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 (A u , B u ), and a nonlinear map Y → G u (Y ), consistent with the nature of the optimisation problem (55).

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 controloriented 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 nonuniform 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.  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 opensource BEM solver Nemoh [5].
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 ∈ [2, 4] [s] (which fully covers the main operational range for this specific device), and a fixed peak-enhancement parameter of γ = 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. η(t), is set to 1600 [s]. Note that such a value corresponds with more than 400 typical wave periods for each possible T w ∈ [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 databased 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.

Data-based control-oriented modelling
Recall that the control-oriented modelling framework, proposed in Sect. 3, is based exclusively on inputoutput data, generated by a suitable numerical nonlinear FK solver. In this paper, we consider the opensource 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 P = {φ j } characterisingf st FK in (13), we choose the trial set which corresponds with the first 16 terms of the polynomial series expansion of f st FK about (z, η) = (0, 0). With respect to the test free-surface elevation η, 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 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 Σ st 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 T w k , 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 P 0 ∈ R 16 , used as a starting point for Algorithm 1, is chosen randomly, Note that the obtained solution exactly coincides with that derived analytically in [36] aiming to highlight the convergence capabilities of the algorithm (see also Remark 11). 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 mappingf st 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 P, explicitly shown in Figure  5 (right), after ≈ 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 η (right, top), and resulting displacement z (right, bottom), together with 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 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 seastates, with H w = 2.5 [m] and T w ∈ [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 Regarding the specific parameters characterising Algorithm 2, the order (dimension) is set toñ = 6, which, as detailed in the reminder of this paragraph, gives a good compromise between computational complexity and model accuracy. The set of (uniformlyspaced) trial time-shifts T c , used to compute an estimate of the (output) time-advance characterising the The results of applying the proposed algorithm are summarised in Figure 9, which presents a Bode plot including each individual ETFE H i ( jω), average ETFEH ( jω) (dotted grey), and frequency-response mapping associated with the computed approximating modelΣ dyn FK (solid black), which clearly presents a good fit with respect to the targetH ( jω). Furthermore, aiming to highlight the difference between standard linear hydrodynamic FK representations, and the proposed approach, Figure 9 also includes the frequencyresponse 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 FK , computed via the nonlinear FK solver Nlkf4all, and that obtained via the control-oriented approximation framework proposed in this study, i.e.ỹ 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 FK −ỹ FK )/ max(|y FK |), consistently showing a satisfactory approximation performance.

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 η for a time window of length T (where energy absorption is to be maximised-see the OCP (37)), we apply a windowing 10 (apodisation) mapping so as to smoothly bring η to zero at the boundaries, and hence the derivative of its periodic extension 10 The particular window mapping considered in this paper is the so-called Planck-taper function [51], which optimally preserves the power spectrum characterising the free-surface elevation (see also the discussion provided in [23]).
is sufficiently smooth. The 'windowed' η 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 finitedimensional 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 t rh seconds, and the procedure is repeated accordingly.

Controller tuning and implementation
The time-window (horizon) length T directly defines the fundamental frequency ω 0 = 2π/T characterising the extended signal generator (39), which, ultimately, defines the approximation space for the optimal con-trol solution f opt u . In particular, the larger the value of T , the smaller the value of ω 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 ι in (55), i.e. the number of harmonics of ω 0 considered to construct the momentbased representation of the optimal energy-maximising control law and to cover the appropriate dynamic frequency range. A large value of ι 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 R ι . In practical scenarios, both T and ι can be tuned together, in terms of a single parameter, i.e. the so-called cutoff frequency ω c = ιω 0 , defining the largest multiple of ω 0 used to construct the extended signal generator (39). In particular, ω 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 ω c fixed, and letting ι = ceil(T ω c /2π), 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 ι = 20, i.e. we consider a total of 10 harmonics of the fundamental frequency ω 0 = 2π/15 [rad/s] in the signal generator (39), and hence the corresponding NP in (55) is carried out over R 20 . Note that this implies a cut-off frequency for the computation of the control solution of ω c ≈ 4 [rad/s].
The receding time step is set to t 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 T 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 t rh , i.e. 0.1 [s]. For the case study presented in this section, the cardinality of T C is hence N 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 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 ≈ 10 −2 [s], which is one order of magnitude smaller than t rh (solid-red line), hence always consistently achieving real-time performance.

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 PI achieving impedance-matching (see e.g. [19,74]) at a (suitably selected) single interpolation frequency. To be precise, K PI : C → C, s → K PI (s), is defined as  Fig. 13 Frequency-response mapping for the optimal control impedance (solid black), and each corresponding K PI (dotted), designed to interpolate I u at each considered peak period T w ∈ [2,4] where the set {θ 1 , θ 2 } ⊂ R can be uniquely selected to interpolate the so-called optimal control impedance I u , at a given interpolation frequency ω I (the reader is referred to "Appendix 3" for the explicit definition of the control impedance I u and the corresponding interpolation procedure).
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 {θ 1 , θ 2 } as a function of the specific operating condition, instead of simply fixing a single controller K PI for all the considered sea-states. In particular, we design (59) to interpolate the optimal control impedance (68) at ω I = 2π/T w , for each considered T w ∈ [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 PI (dotted), designed to interpolate I u for each corresponding (input) seastate condition. Note that, in (natural) resonance conditions (T w ≈ 3 [s]), the optimal K PI is simply a constant value, i.e. it becomes passive.  Figure 14 displays the main performance results obtained with the proposed nonlinear optimal momentbased 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).
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.

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 (bidirectional) 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.
Acknowledgements The results of this publication reflect only the author's view and the European Commission is not responsible for any use that may be made of the information it contains Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No 101024372.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

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 for any A 1 ∈ R k×l , A 2 ∈ R l×m , and A 3 ∈ R m×n , holds. An application of the vec operator to Eq. (48), together with the property defined via (60) where we have used the fact that S = −S . To finalise proving the claim, we note that the map (Π, L η ) → G(Π, L η ) can be fully characterised in terms of Y , i.e. the map 11 (Y , L η ) → G(Y , L η ), by virtue of the following relation. Recall that the state-vector x of system (31) is x = [z,ż] = [z y], and that, in steady-state conditions, x ss = Π ξ and y ss = CżΠξ = Y ξ . Following [25,69], one can derive that and hence Π can be fully written in terms of Y as 11 Though the domain of the map actually changes from R 2×ι × R 1×ι to R 1×ι × R 1×ι , we keep the notation G for the sake of clarity and convenience.

Appendix 2: On the derivation of Eq. (56)
Recall the definition of the state and input constraints C in (36). In steady-state conditions, one can explicitly use the approximated moment Y to map the set C as C : where we have explicitly used the relation posed in equation (63). Using the property that, for any f ∈ C 1 and α ∈ R + /0, | f (t)| ≤ α → f (t) ∈ [−α, α], and enforcing the set C ss at the collocation set T C ⊂ R + , we can write (65) as from which the pairs of matrices (A z , B z ) and (Aż, Bż) in (57) follow directly. Finally, the expressions for the pair (A u , B u ), and the nonlinear map Y → G u (Y ) in (57), can be obtained by noting that for any fixed L η , via Remark 24.