# Spectral Control of Wave Energy Converters with Non-Ideal Power Take-off Systems

Alexis Mérigaud, Paolino Tona

► To cite this version:

Alexis Mérigaud, Paolino Tona. Spectral Control of Wave Energy Converters with Non-Ideal Power Take-off Systems. Journal of Marine Science and Engineering, 2020, 8 (11), pp.851. 10.3390/jmse8110851. [https://doi.org/10.3390/jmse8110851](https://doi.org/10.3390/jmse8110851)

HAL Id: hal-03097665

https://ifp.hal.science/hal-03097665

Submitted on 5 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Spectral Control of Wave Energy Converters with Non-Ideal Power Take-off Systems

Alexis Mérigaud * and Paolino Tona

IFP Energies Nouvelles, 1-4 Avenue du Bois Préau, 92852 Rueil-Malmaison, France; paolino.tona@ifpen.fr

* Correspondence: alexis.mérigaud@ifpen.fr

Received: 7 October 2020; Accepted: 21 October 2020; Published: 28 October 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Spectral control is an accurate and computationally efficient approach to power-maximising control of wave energy converters (WECs). This work investigates spectral control calculations with explicit derivative computation, applied to WECs with non-ideal power take-off (PTO) systems characterised by an efficiency factor smaller than unity. To ensure the computational efficiency of the spectral control approach, it is proposed in this work to approximate the discontinuous efficiency function by means of a smooth function. A non-ideal efficiency function implies that the cost function is non-quadratic, which requires a slight generalisation of the derivative-based spectral control approach, initially introduced for quadratic cost functions. This generalisation is derived here in some detail given its practical interest. Two application case studies are considered: the Wavestar scale model, employed for the WEC control competition (WECCCOMP), and the 3rd reference model (RM3) two-body heaving point absorber. In both cases, with the approximate efficiency function, the spectral approach calculates WEC trajectory and control force solutions, for which the mean electrical power is shown to lie within a few percent of the true optimal electrical power. Regarding the effect of a non-ideal PTO efficiency upon achievable power production, and concerning heaving point-absorbers, the results obtained are significantly less pessimistic than those of previous studies: the power achieved lies within 80–95% of that obtained by simply applying the efficiency factor to the optimal power with ideal PTO.

**Keywords:** wave energy converters; spectral control; optimal control; power take-off

## 1. Introduction

Power-maximising control is a promising path to improve the economic competitiveness of WECs [1]. Various control methods have been proposed, of which the most basic are simple parametric control strategies, whereby one or two scalar control parameters are adjusted on a sea state-by-sea state basis. Recent years have witnessed a growing interest for model predictive control (MPC) strategies [2], whereby an optimisation problem is solved in real-time to update a prescribed control input, or a reference trajectory to be followed by the WEC. MPC strategies take profit from wave excitation predictions. They can, in theory, handle nonlinear WEC models, through the use of nonlinear optimisation algorithms, as well as constraints on the WEC dynamics, control force or instantaneous power.

First introduced in [3] for WEC control applications, spectral and pseudo-spectral MPC methods discretise the wave inputs and optimisation variables by means of basis functions, typically Fourier [4,5] (when a periodic wave input is considered), or Chebyshev-like [6] (to handle directly the non-periodic wave signal, seen within the receding-horizon control window). In terms of signal approximation (i.e., when it comes to approximate a given function by means of a finite number of “spectral” basis functions), spectral techniques are said to have *superlinear convergence* properties [7], which means that the approximation error decreases more than linearly with the number of basis functions employed,
provided that the function approximated is continuous (with a geometric convergence rate for functions with a finite degree of smoothness, and an exponential convergence rate for infinitely-smooth functions such as linear functions). As a consequence, spectral methods require a small number of coefficients to describe the control inputs or WEC dynamics with a prescribed accuracy, thus reducing the number of optimisation variables, and the associated computation cost. Spectral methods can thus handle relatively long control horizons. In addition, the “aggressiveness” of the control solution, i.e., how abruptly the solution oscillates, can be tuned, by adjusting the number of basis functions.

In [5], a computationally efficient spectral control formulation for nonlinear, 1-DoF WEC models is introduced in which, using the dynamical equation, the control variable is eliminated in the expression of the objective function, thus resulting in a problem where the only variables to optimise are those, corresponding to the WEC motion. The approach in [5], combined with the explicit computation of the first- and second-order derivatives of the objective function, is shown to be orders of magnitude faster than usual spectral control formulations, whereby the dynamical equations are considered as a set of equality constraints. Although the proposed approach is based on Fourier basis functions, it has been successfully implemented in a receding-horizon set-up [8,9]. Also note that combining spectral control (for optimal reference generation) with tracking control entails a certain degree of robustness to modelling errors, in contrast to direct computation of the optimal control force, which requires accurate modelling of all the WEC dynamics. Indeed, the controller presented in [9] is shown to be little affected by tracking loop imperfection as well as by estimation and forecast errors.

In all references previously mentioned, like in the vast majority of the wave energy control literature, only hydrodynamic power absorption is considered. From an optimal control perspective, the objective function, here denoted as $f$, is generally the raw mechanical power absorbed by the PTO system [2]. For a WEC with one mechanical degree of freedom, the raw absorbed power is written as follows,

$$f(x, u) = -\frac{1}{T} \int_{t=0}^{T} \dot{x}udt \qquad (1)$$

where $T$ is the length of the control horizon, $x$ denotes the WEC generalised position and $u$ the PTO force (which is the control input).

However, in practical applications, the PTO system is not ideal, so that the electrical power, $P_e$, is not equal to the absorbed mechanical power, $P_a = -u\dot{x}$, as highlighted by a number of authors [10–17]. More specifically, when power is mechanically absorbed from the waves, the electrical power $P_e$, actually available at the end of the mechanical-to-electrical conversion stage, is smaller than the absorbed power $P_a$, because of the losses occurring in the conversion stage ($0 \le P_e \le P_a$). Conversely, if a reactive control strategy is adopted [1], at times some power must be returned from the electromechanical system back into the ocean ($P_a \le 0$). At those instants, because of the losses in the PTO system, the electrical power provided by the grid to the PTO system must be greater than $|P_a|$, so that $P_e \le P_a \le 0$. PTO losses should be taken into account when designing the WEC control strategy [10–17]. In fact, a reactive control strategy, obtained without modelling losses in the objective function, can yield negative average generated power when applied to the true system, where losses occur (see, for example, in [13]).

Therefore, it is of high practical interest to extend the methodology proposed in [5] to more general objective functions of the form

$$f(x, u) = \frac{1}{T} \int_{t=0}^{T} \alpha(x, \dot{x}, u)dt \qquad (2)$$

where $\alpha$ is some scalar function of $x, \dot{x}$ and $u$ (while the control calculations presented in [5] can only deal with cost functions in the form $-u\dot{x}$, quadratic in the problem variables). $P_e$, being modelled as a function of $P_a = -u\dot{x}$, is a special case of the formulation in Equation (2). Such a formulation can also include, for example, quadratic terms introduced to penalise the WEC excursions or control effort.
This work, like that in [5], is concerned with the calculation of optimal control solutions in an “off-line” fashion (as opposed to receding-horizon implementations such as in [6,8,9]), taking into account the totality of a relatively long, simulated irregular wave input signal comprising a few tens of wave periods. Therefore, the insight provided, regarding the effect of a non-ideal PTO system onto optimal control results, is reasonably general, and remains free of the many effects, which interact when considering a particular receding-horizon implementation [9]. Furthermore, the optimal control solution, obtained off-line, and the corresponding generated power, provide a reference to assess the performance of real-time, receding-horizon control algorithms (which cannot, in theory, outperform the optimal solution obtained off-line). The present study is restricted to WECs with one degree of freedom.

The rest of this paper is organised as follows.

* Section 2 revisits the method presented in [5], for an objective function of the form of Equation (2).

* Section 3 introduces the simple non-ideal PTO model retained in this work, and two numerical case studies: the Wavestar model used in the WECCCOMP challenge [18], and the RM3 two-body heaving point absorber [19].

* For both case studies, numerical results are shown and discussed in Section 4;

* Finally, conclusions are presented in Section 5.

## 2. An Efficient Spectral Control Method for a 1-DoF, Non-Linear WEC Model

### 2.1. WEC Dynamical Model

Consider a WEC, modelled with one degree of freedom, denoted $x$. Newton’s second law, which describes the WEC motion, typically takes the following mathematical form (see in [5], and detailed justifications in Chapter 4 of the work in [20]):

$$ \mathcal{L}\{x\} + n(x, \dot{x}) = u + e \quad (3) $$

where the dependence in the time variable $t$ has been omitted to alleviate notations, and

* $\mathcal{L}$ represents a linear integro-differential operator, containing some of, or all the terms of the well-known Cummins’ equation [21]:

$$ \mathcal{L}\{x\}(t) := (I + I_{\infty})\ddot{x}(t) + C\dot{x}(t) + \int_{\tau=0}^{\infty} k_r(\tau)\dot{x}(t-\tau)d\tau + Sx(t) \quad (4) $$

where $I$ represents the system inertia; $I_{\infty}$ and $k_r$ are the radiation infinite-frequency added inertia and damping convolution kernel, respectively; C models some linear damping term (for example, to take friction effects into account); and $S$ is the hydrostatic restoring force (resulting from the balance between gravity and Archimede’s force);

* $n$ represents some nonlinear function, with the aim of modelling nonlinear dynamical effects, not accounted for by a linear modelling approach. If $n$ includes effects, usually modelled linearly in Cummins’ Equation (4), then the corresponding linear coefficients in Equation (4) can be assumed equal to zero;

* $e(t)$ is the wave excitation force; and

* $u(t)$ is the control force, exerted by the PTO system onto the WEC’s moving part.

The operator $\mathcal{L}$ is represented by a frequency-response mapping $z(\omega) = S - \omega^2(I + A_r(\omega)) - j\omega(C + B_r(\omega))$, where the radiation frequency-dependent added mass $A_r(\omega)$ and damping $B_r(\omega)$ are related to $I_{\infty}$ and $k_r$ through Ogilvie’s relation [22].
## 2.2. Spectral Control Problem Formulation

Consider an optimal control problem of the following form,

$$ \begin{cases} \max_{x,u} f(x, u) \\ \text{s.t. } \mathcal{L}\{x\} + n(x, \dot{x}) = u + e \end{cases} \qquad (5) $$

where $f(x, u) := \frac{1}{T} \int_{t=0}^{T} \alpha(x, \dot{x}, u)dt$; T is the control time horizon; the equality constraint expresses the fact that the dynamical Equation (3) must be satisfied at every instant; and the function $\alpha$, as mentioned in the introduction, represents the “instantaneous” objective function. Given the particular form of the dynamical equation, u can be expressed as a function of the other variables, so that Problem (5) can be recast as follows,

$$ \max_x f(x) \qquad (6) $$

where $f(x) = \frac{1}{T} \int_{t=0}^{T} \alpha(x, \dot{x}, \mathcal{L}\{x\} + n(x, \dot{x}) - e)dt$, while the variable u, as well as the equality constraint, are eliminated, so that only the device trajectory is now optimised.

Assume a periodic, polychromatic excitation signal with period T, of the form

$$ e(t) = \frac{1}{\sqrt{2}} \hat{e}_1 + \sum_{n=1}^{N} \hat{e}_{2n} \cos(\omega_n t) + \hat{e}_{2n+1} \sin(\omega_n t) \qquad (7) $$

where the frequencies $\omega_n$ are defined harmonically with $\Delta\omega = \frac{2\pi}{T}$ and $\forall n \in \{1...N\}, \omega_n = n\Delta\omega$. In the following, consistently with a periodic excitation signal, solutions x are searched among periodic signals with period T, and the WEC motion is assumed of the form

$$ x(t) = \frac{1}{\sqrt{2}} \hat{x}_1 + \sum_{n=1}^{N} \hat{x}_{2n} \cos(\omega_n t) + \hat{x}_{2n+1} \sin(\omega_n t) \qquad (8) $$

Note that, in Equations (7) and (8), the number N of harmonics being formally identical does not imply any loss of generality: for example, the number of harmonics in x can be considered larger than that in e, by setting to zero higher frequency components of e. In the extreme case of a zero-mean, monochromatic wave input, all components $\hat{e}_n$, except for $\hat{e}_2$ and $\hat{e}_3$, are zero, while the solution x can have an arbitrarily large number of components N.

The next step consists of discretising the integral in objective function, using a set of M equally-spaced points $t_m$, spanning the interval $[0; T]$, with typically $M \ge 2N + 1$, so that the maximisation problem becomes

$$ \max \tilde{f}(x) := \sum_{m=1}^{M} \alpha(t_m) \qquad (9) $$

where, for the sake of conciseness, $\alpha(t_m)$ denotes the value of $\alpha$ when its arguments are evaluated at time $t_m$.

Define the matrix $\Phi \in \mathbb{R}^{M \times (2N+1)}$ as

$$ \forall i \in \{1...M\}, j \in \{1...N\}, \begin{cases} \Phi_{i,1} = \frac{1}{\sqrt{2}} \\ \Phi_{i,2j} = \cos(\omega_j t_i) \\ \Phi_{i,2j+1} = \sin(\omega_j t_i) \end{cases} \qquad (10) $$

and the matrix

$$ \Omega = \begin{pmatrix} \Omega_0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \Omega_N \end{pmatrix} \qquad (11) $$
where $\Omega_0 = 0$ and, $\forall n \ge 1$, $\Omega_n = \begin{pmatrix} 0 & \omega_n \\ -\omega_n & 0 \end{pmatrix}$.

Define vectors $\hat{x} := (\hat{x}_1 \cdots \hat{x}_N)^T$, $\hat{e} := (\hat{e}_1 \cdots \hat{e}_N)^T$, $\mathbf{x} := (x(t_1) \cdots x(t_M))^T$, $\dot{\mathbf{x}} := (\dot{x}(t_1) \cdots \dot{x}(t_M))^T$ and $\mathbf{e} := (e(t_1) \cdots e(t_M))^T$. Then, it is easy to see that $\mathbf{x} = \Phi\hat{x}$, $\dot{\mathbf{x}} = \Phi\Omega\hat{x}$ and $\mathbf{e} = \Phi\hat{e}$.

Furthermore, the linear terms of $\mathcal{L}$ evaluated at times $t_1...t_M$, $l := (\mathcal{L}\{x\}(t_1) \cdots \mathcal{L}\{x\}(t_M))^T$, can be calculated as $l = \mathbf{L}\hat{\mathbf{x}}$, where

$$ \mathbf{L} = \begin{pmatrix} L_0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & L_N \end{pmatrix} \qquad (12) $$

with blocks $L_0 = \Re\{z(0)\}$ and $\forall n \in \{1...N\}$, $L_n = \begin{pmatrix} \Re\{z(\omega_n)\} & \Im\{z(\omega_n)\} \\ -\Im\{z(\omega_n)\} & \Re\{z(\omega_n)\} \end{pmatrix}$.

Finally, the objective function can be expressed as a function of $\hat{x}$, as follows,

$$ \tilde{f}(\hat{x}) = \mathbf{1}_{M \times 1} \alpha(\Phi\hat{x}, \Phi\Omega\hat{x}, \Phi\mathbf{L}\hat{x} + n(\Phi\hat{x}, \Phi\Omega\hat{x}) - \Phi\hat{e}) \qquad (13) $$

where $\mathbf{1}_{M \times 1}$ is the unit vector of size $M \times 1$ and, with a slight abuse of notation, the functions $\alpha$ and $n$ are applied component-wise to the vectors used as arguments, i.e., for example, $g(\mathbf{v}) = (g(\mathbf{v}_1)...g(\mathbf{v}_M))^T$.

Note that Equation (13) differs from the objective function in [5], as $\alpha$ is, in general, not a quadratic function of its arguments. The calculations which follow constitute a useful generalisation of the results in [5].

## 2.3. Explicit Calculation of the Objective Function Derivatives

In this section, it is assumed that $\alpha$ and $n$ are continuous in all their arguments, and admit second-order derivatives. Under such an assumption, the first- and second-order derivatives of $\tilde{f}$, with respect to the components of $\hat{x}$, can be put to good use within a gradient-based optimisation algorithm [5]. In the following, it is shown how to calculate explicitly those derivatives. The calculations do not present any particular difficulty, but they are not detailed here for the sake of brevity.

A useful convention is first introduced. Given a function $y$ defined on $[0; T]$, the notation $D_y$ represents the $M \times M$ diagonal matrix with diagonal entries $y(t_1), ...y(t_M)$. Then, the vector of 1st-order derivatives is given as

$$ \frac{\partial \tilde{f}}{\partial \hat{\mathbf{x}}} = \mathbf{1}_{M \times 1} [D_{v_1} \Phi + D_{v_2} \Phi \Omega + D_{v_3} \Phi \mathbf{L}] \qquad (14) $$

with

$$ \left\{ \begin{array}{l} v_1 = \frac{\partial \alpha}{\partial x} + \frac{\partial \alpha}{\partial u} \frac{\partial n}{\partial x} \\ v_2 = \frac{\partial \alpha}{\partial \hat{x}} + \frac{\partial \alpha}{\partial u} \frac{\partial n}{\partial \hat{x}} \\ v_3 = \frac{\partial \alpha}{\partial u} \end{array} \right. $$

The matrix of second-order derivatives is given as

$$ \begin{aligned} \frac{\partial^2 \tilde{f}}{\partial \hat{x}^2} = & \Phi^\top D_{w_1} \Phi + \Phi^\top D_{w_2} \Phi \Omega + \Omega^\top \Phi^\top D_{w_2} \Phi \\ & + \Omega^\top \Phi^\top D_{w_3} \Phi \Omega + \Phi^\top D_{w_4} \Phi \mathbf{L} + \mathbf{L}^\top \Phi^\top D_{w_4} \Phi \\ & + \Omega^\top \Phi^\top D_{w_5} \Phi \mathbf{L} + \mathbf{L}^\top \Phi^\top D_{w_5} \Phi \Omega \\ & + \mathbf{L}^\top \Phi^\top D_{w_6} \Phi \mathbf{L} \end{aligned} \qquad (15) $$
with:

$$
\left\{
\begin{aligned}
w_1 &= \frac{\partial^2 \alpha}{\partial x^2} + 2 \frac{\partial^2 \alpha}{\partial x \partial u} \frac{\partial n}{\partial x} + \frac{\partial^2 \alpha}{\partial u^2} \left( \frac{\partial n}{\partial x} \right)^2 + \frac{\partial \alpha}{\partial u} \frac{\partial^2 n}{\partial x^2} \\
w_2 &= \frac{\partial^2 \alpha}{\partial \dot{x} \partial x} + \frac{\partial^2 \alpha}{\partial u^2} \frac{\partial n}{\partial x} \frac{\partial n}{\partial \dot{x}} + \frac{\partial^2 \alpha}{\partial \dot{x} \partial u} \frac{\partial n}{\partial x} + \frac{\partial^2 \alpha}{\partial x \partial u} \frac{\partial n}{\partial \dot{x}} + \frac{\partial \alpha}{\partial u} \frac{\partial^2 n}{\partial x \partial \dot{x}} \\
w_3 &= \frac{\partial^2 \alpha}{\partial \dot{x}^2} + 2 \frac{\partial^2 \alpha}{\partial \dot{x} \partial u} \frac{\partial n}{\partial \dot{x}} + \frac{\partial^2 \alpha}{\partial u^2} \left( \frac{\partial n}{\partial \dot{x}} \right)^2 + \frac{\partial \alpha}{\partial u} \frac{\partial^2 n}{\partial \dot{x}^2} \\
w_4 &= \frac{\partial^2 \alpha}{\partial u \partial x} + \frac{\partial^2 \alpha}{\partial u^2} \frac{\partial n}{\partial x} \\
w_5 &= \frac{\partial^2 \alpha}{\partial u \partial \dot{x}} + \frac{\partial^2 \alpha}{\partial u^2} \frac{\partial n}{\partial \dot{x}} \\
w_6 &= \frac{\partial^2 \alpha}{\partial u^2}
\end{aligned}
\right.
$$

In practice, at each iteration of the optimisation algorithm, the objective function $\tilde{a}$, the nonlinear terms $n$, and their 1st- and 2nd-order derivatives, are evaluated along the trajectory given by $\mathbf{x} = \Phi\dot{\mathbf{x}}$. $\tilde{f}$ is evaluated; $v_1...v_3$ are built in order to calculate $\frac{\partial \tilde{f}}{\partial \dot{x}}$ as in Equation (14), and $w_1...w_6$ are built in order to calculate $\frac{\partial^2 \tilde{f}}{\partial \dot{x}^2}$ as in Equation (15).

## 2.4. Inequality Constraints

The reader may have noticed that inequality constraints have not been introduced in the problem formulation. Yet, operational limitations, typically on the WEC position or velocity, on the control force or on instantaneous power, are essential for the safe operation of the device.

In particular, constraints which are nonlinear with respect to $x$ (e.g., force or power constraints) are computationally challenging. How their 1st- and 2nd-order derivatives can be explicitly computed, is detailed in [5]. The more general objective function, introduced in the present paper in Equation (5), does not result in any change in the constraint-related calculations, which are therefore not reproduced here.

# 3. Numerical Case Studies

## 3.1. Non-Ideal Power Take-Off Systems

The effect of a non-ideal PTO efficiency upon WEC control is specifically addressed in a number of studies [10–15,18]. In most of them [11,12,14–18], the non-ideal PTO system is modelled by means of a non-unity efficiency factor $\mu$ verifying $0 \le \mu \le 1$, so that the instantaneous electric power is expressed as follows.

$$P_e = \begin{cases} -\mu u \dot{x} & \text{if } -u \dot{x} \ge 0 \\ -\frac{1}{\mu} u \dot{x} & \text{if } -u \dot{x} < 0 \end{cases} \quad (16)$$

The simple model (16) is also retained in the present study, for the two WEC models described in Sections 3.2 and 3.3.

Reactive control techniques (which imply that, at times, $P_a$ is negative) are essential to achieve optimal or near-optimal hydrodynamic power absorption, when the incoming wave frequency is away from the system's resonant frequency. In the case of monochromatic waves, a simple PI control law, of the form $u = -b\dot{x} - kx$, can achieve optimal power absorption, following the principle of impedance matching [1]. In irregular waves, a PI control law also presents significant improvements in absorbed power, with respect to a passive law of the form $u = -b\dot{x}$, although the PI controller is then sub-optimal. A salient point of studies [11,12] is that reactive control is particularly affected by a non-ideal PTO efficiency. In particular, it appears that the smaller $\mu$ is, the less reactive control should be performed (i.e., $|k|$ should be smaller).
Overall, all studies [10–15,17,18] are relatively pessimistic, regarding the mean electrical power $\bar{P}_e$ which can be achieved with a non-ideal PTO system: a 10% decrease in efficiency results in a drop in $\bar{P}_e$ of more than 50%, even when control parameters are tuned taking the PTO efficiency into account<fnref n="1" />.

Define $h(P_a)$ the efficiency function, so that $h(P_a) = \frac{P_e}{\bar{P}_a}$. Consistently with Equation (16), $h$ is equal to $\mu$ for positive $P_a$, and $1/\mu$ for negative $P_a$, as illustrated in Figure 1. Such a discontinuous function is undesirable for the efficient implementation of the spectral control approach, detailed in Section 2, and which requires the explicit calculation of the 1st- and 2nd-order derivatives of the objective function. Therefore, in this work, $h(P_a)$ is approximated by the following function,

$$h_{\kappa}(P_a) = A \tanh(\kappa P_a) + B \qquad (17)$$

where $A = \frac{1}{2}(\mu - 1/\mu)$, $B = \frac{1}{2}(\mu + 1/\mu)$ and $\kappa > 0$ is a real parameter, which governs the accuracy of the approximation, as illustrated in Figure 1. Note that in [16], the use of a smoothed approximation of the function (16) is also advocated, in order to avoid issues in gradient-based optimisation for the considered optimal control problems. However, the sensitivity of the results to smoothing has not been investigated.

<figure><img src="image_4.png" /><figcaption>Figure 1. Efficiency function $h$ for $\mu = 0.9$ and $\mu = 0.7$, and their approximations $h_{\kappa}$.</figcaption></figure>

If the exact efficiency is considered, the function $\alpha$ in the optimisation problem (5) is $\alpha(\dot{x}, u) := -u\dot{x}h(-u\dot{x})$. However, in practice, the problem is solved using the approximate function $\alpha_{\kappa}(\dot{x}, u) := -u\dot{x}h_{\kappa}(-u\dot{x})$. For all $(\dot{x}, u)$, $\alpha(\dot{x}, u) \le \alpha_{\kappa}(\dot{x}, u)$. Therefore, the optimal mean electric power, $\bar{P}_{e,\kappa}^*$, obtained with $\alpha_{\kappa}$, is an upper bound to the optimal exact electric power, $\bar{P}_{e,\kappa}^*$, which would be obtained if Problem (1) could be solved using the exact efficiency. Furthermore, for all $\kappa_1 \le \kappa_2$, for all $(\dot{x}, u)$, $\alpha_{\kappa_2}(\dot{x}, u) \le \alpha_{\kappa_1}(\dot{x}, u)$, which implies that $\bar{P}_{e,\kappa}^*$ decreases monotonically with $\kappa$, and has $\bar{P}_{e}^*$ as a limit when $\kappa \to \infty$.

$\bar{P}_{e,\kappa}^*$ can also be bounded from below: consider the solution $x_{\kappa}^*, u_{\kappa}^*$, obtained with the approximate objective function. $x_{\kappa}^*, u_{\kappa}^*$ is necessarily sub-optimal, with respect to the exact objective function. Therefore, the exact mean electric power, effectively achieved with $(x_{\kappa}^*, u_{\kappa}^*)$, constitutes a lower bound to $\bar{P}_{e,\kappa}^*$, and is calculated as follows.

$$\bar{P}_{e,\kappa}^{\dagger} = \frac{1}{T} \int_{t=0}^{T} \alpha(\dot{x}_{\kappa}^{*}, u_{\kappa}^{*}) dt \qquad (18)$$

<footnote n="1">¹ It should be noted, however, that the control strategies adopted, and the PTO efficiency models, differ across studies, so that accurate comparisons are difficult.</footnote>
Therefore the following inequality holds.

$$ \bar{P}_{e,\kappa}^{\dagger} \leq \bar{P}_{e}^{*} \leq \bar{P}_{e,\kappa}^{*} \qquad (19) $$

Even if the true optimal electrical power $\bar{P}_{e}^{*}$ cannot be calculated for the exact efficiency, by setting $\kappa$ to a large enough value, it can be ensured that the electric power $\bar{P}_{e,\kappa}^{\dagger}$ effectively achieved with the solution found, lies within a prescribed percentage of $\bar{P}_{e}^{*}$.

Finally, let $(x^*, u^*)|_{\mu}$ be the optimal solution for efficiency $\mu < 1$, and $\bar{P}_{e}^{*}|_{\mu}$ (resp. $\bar{P}_{a}|_{\mu}$) the mean electrical (resp. absorbed) power with the solution $(x^*, u^*)|_{\mu}$. It is easy to show<fnref n="2" /> that $\bar{P}_{e}^{*}|_{\mu} < \mu \bar{P}_{a}|_{\mu}$. Define $\bar{P}_{id}^{*}$ as the optimal power for an ideal PTO (absorbed and electric power are then identical). The solution $(x^*, u^*)|_{\mu}$ is sub-optimal with respect to an ideal PTO, thus $\bar{P}_{a}|_{\mu} < \bar{P}_{id}^{*}$. Therefore, the following simple inequality holds.

$$ \bar{P}_{e}^{*}|_{\mu} < \mu \bar{P}_{id}^{*} \qquad (20) $$

Equation (20) means that the optimal mean electric power, for a non-ideal PTO efficiency, cannot exceed the optimal mean electric power for an ideal PTO efficiency, multiplied by the efficiency factor. In the particular case, where the optimal solution, for an ideal PTO efficiency, is purely passive, then that solution is also the optimal solution for a non-ideal PTO efficiency $\mu < 1$, and the inequality becomes an equality: $P_{e}^{*}|_{\mu} = \mu P_{e}^{*}|_{id}$. This can happen if the incoming wave frequency coincides with the WEC natural frequency, so that no reactive power is needed to achieve optimal control.

## 3.2. The WECCCOMP Model

The Wavestar device is a point absorber WEC, primarily operating in heave, rigidly connected to a rotation arm. The power is captured through the arm rotation around a fixed axis. In the WEC control competition (WECCCOMP), a scale model of the Wavestar device is considered. The linearised WEC dynamics can be modelled by expressing Newton's second law around the rotation axis, which results in a dynamical equation, of the form of Equation (3), where the nonlinear terms $n$ are set to zero. The variable $x$, in this case, corresponds to the angular displacement of the rotation arm. All the linear parameters of Equation (3) are derived from the linear model, provided by the competition organisers [18]. Lower and upper limits on the PTO torque (±10 N.m) and angular position (±24°) are considered as constraints. Six JONSWAP wave spectra [23], proposed by the organisers [18], and with characteristics listed in Table 1, are used to generate random excitation torque time series, with duration $T = 40$ s, at least equal to 20 typical wave periods for the sea states considered.

In each random realisation, the optimal solution $\hat{x}_{\kappa}^{*}|_{\mu}$, and power $\bar{P}_{e,\kappa}^{*}|_{\mu}$, with $\mu = 1, 0.9, \dots, 0.5$, are calculated using the *approximate* efficiency function $h_{\kappa}$. Then, the *exact* power $\bar{P}_{e,\kappa}^{\dagger}|_{\mu}$, obtained with $\hat{x}_{\kappa}^{*}|_{\mu}$, is also calculated.

For each pair of sea states (1–4, 2–5, 3–6), the appropriate range of values for $\kappa$ differs (when typical values for absorbed power are larger, a smaller $\kappa$ yields a satisfactory approximation). For each pair of sea states, $\kappa$ is gradually increased, until bounds given by Equation (19) are sufficiently close. Suitable values for $\kappa$ are indicated in Table 1.

<table><caption>**Table 1.** List of WECCCOMP JONSWAP sea states.</caption><thead><tr><th>ID</th><th>(H<sub>m0</sub>, T<sub>p</sub>, γ)</th><th>ID</th><th>(H<sub>m0</sub>, T<sub>p</sub>, γ)</th><th>κ</th></tr></thead><tbody><tr><td>1</td><td>(0.0208, 0.988, 1)</td><td>4</td><td>(0.0208, 0.988, 3.3)</td><td>100</td></tr><tr><td>2</td><td>(0.0625, 1.412, 1)</td><td>5</td><td>(0.0625, 1.412, 3.3)</td><td>5</td></tr><tr><td>3</td><td>(0.1042, 1.836, 1)</td><td>6</td><td>(0.1042, 1.836, 3.3)</td><td>2</td></tr></tbody></table>

<footnote n="2">² To see this, decompose $\bar{P}_{a}$ into its passive and reactive components.</footnote>
### 3.3. The RM3 Device

The RM3 device is a self-reacting, two-body point absorber, exploiting the relative heave motion of the two bodies [19]. The RM3 consists of a cylindrical float (B1) which vertically oscillates with waves, relative to a vertical column spar buoy (B2) damped by a subsurface reaction plate. Within the scope of linear hydrodynamics, the heaving mode of motion is decoupled from the other degrees of freedom, and therefore B1 and B2 may be modelled in terms of heave only. Assuming harmonic signals and adopting complex notations (omitting the frequency dependence for conciseness), the position in heave of B1 and B2, $x_1$ and $x_2$, can be described by the following two equations [24],

$$ \begin{cases} z_1\hat{x}_1 + z_{21}\hat{x}_2 = \hat{e}_1 + \hat{u} \\ z_{12}\hat{x}_1 + z_{22}\hat{x}_2 = \hat{e}_2 - \hat{u} \end{cases} \qquad (21) $$

where

* the complex mechanical impedance $z_i$ is built as $z_i(\omega) = S_i - \omega^2(I_i + A_{r,i}(\omega)) - j\omega(C_i + B_{r,i}(\omega))$ from the hydrodynamic coefficients of Body i,

* the coupling terms $z_{ij}(\omega) = z_{ji}(\omega) = -\omega^2 A_{r,i \to j}(\omega) - j\omega B_{r,i \to j}(\omega)$ are due to the radiation interaction between B1 and B2,

* $\hat{e}_i$ denotes the wave excitation force onto Body i and

* the PTO system applies a force $u$ from B2 onto B1.

By rearranging Equation (21), it is easy to obtain a new formulation in terms of the relative motion<fnref n="3" />, $x := x_1 - x_2$, sufficient to solve the control problem (see in [24] for more detail):

$$ z_{eq}\hat{x} = \hat{e}_{eq} + \hat{u} \qquad (22) $$

where $z_{eq} = \frac{z_2z_1-z_{12}^2}{z_0}$, $\hat{e}_{eq} = \frac{z_2+z_{12}}{z_0}\hat{e}_1 + \frac{z_1+z_{12}}{z_0}\hat{e}_2$ and $z_0 = z_1 + z_2 + 2z_{12}$. The proposed spectral method can thus be applied, using the equivalent impedance $z_{eq}$ to determine the entries of the matrix L as in Section 2.2.

Control calculations are carried out in four JONSWAP [23] sea states, inspired by the scatter diagram of the device's design location, given in [19], with $H_{m0} = 2m$, $\gamma = 3$, and $T_p = 6, 8, 10$ and $12$ s. The relative heave $x$ is subject to the constraint $-2m \le x \le 2m$ [19]. The wave signal duration is set to $T = 300$ s (more than 25 typical wave periods in all sea states). A value $\kappa$ of the order of $10^{-6}$ allows for accurate optimal power estimates<fnref n="4" />.

## 4. Numerical Results

All constrained control results were obtained using the interior-point algorithm implemented in the Matlab<fnref n="5" /> `fmincon` function. For the WECCCOMP model, the cut-off frequency for the control solution is set to 4 Hz (as the signal duration is $T = 40$ s, the order $N$ of the approximation is therefore $N = 4 \times 40 = 160$). The number $M$ of collocation points is set to $M = 3N + 1$. For the RM3 device, the cut-off frequency for the control solution is set to 0.5 Hz (as the signal duration is $T = 300$ s, the order $N$ of the approximation is therefore $N = 0.5 \times 300 = 150$). The number $M$ of collocation points is set to $M = 4N + 1$.

<footnote n="3">³ For the motion to be fully determined, another equation is necessary, for example in terms of $y := x_1 + x_2$, but this does not have any influence in terms of power absorption.</footnote>

<footnote n="4">⁴ The reader may have noted that this value is very different from those, which were found appropriate for the WECCCOMP model. This is because, given the form of Equation (17), the appropriate value for $\kappa$ depends on the typical magnitude of $P_a$ in the RM3 case, which is a full-scale model, $P_a$ takes much larger values than in the WECCCOMP case.</footnote>

<footnote n="5">⁵ www.matlab.com.</footnote>
## 4.1. The WECCCOMP Model

Figure 2 shows control results obtained with the WECCCOMP numerical model (see Section 3.2), in the six sea states of the competition, and values of $\mu$ between 1 and 0.5. In each sea state, power values are normalised by $\bar{P}_{id}^*$, the optimal power obtained with $\mu = 1$. $\bar{P}_{id}^*$ essentially represents the optimal absorbed mechanical power; it is a well-known theoretical limit in the wave energy literature [1], which can be calculated as $\bar{P}_{id}^* = \int_{\omega=0}^{\infty} \frac{-\omega|\hat{e}(\omega)|^2}{2\Im\{z(\omega)\}} d\omega$. As per Equation (19), $\bar{P}_{e,\kappa}^\dagger|_{\mu}$ and $\bar{P}_{e,\kappa}^*|_{\mu}$ together provide bounds to $\bar{P}_e^*|_{\mu}$, the optimal power for the exact efficiency function.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Normalised mean electric power for the model Wavestar device in the six WECCCOMP sea states, over a range of efficiency parameter values.</figcaption></figure>

Also shown in Figure 2, the dotted lines labelled "$\mu\bar{P}_{id}^*$" represent the upper bound of Equation (20), i.e., a hypothetical scenario where $\bar{P}_{id}^*$ would "only" be affected by a factor of $\mu$. $\bar{P}_{mism.}$ is the actual electric power which would be obtained, if the solution $\hat{x}_{id}^*$ to the ideal control problem were followed by a WEC with imperfect PTO efficiency $\mu \le 1$; in other words, $\bar{P}_{mism.}$ shows the effect of a *mismatch* between the nominal efficiency ($\mu = 1$) and the actual efficiency ($1 > \mu \ge 0.5$). Finally, $\bar{P}_{Bpto}$ represents the power obtained using a simple linear PTO damping, with coefficient $B_{pto}$ optimised for each sea state.

The importance of taking into account the non-ideal PTO efficiency in control calculations is evidenced by $\bar{P}_{mism.}$, which can even take negative values, similarly to the results of [13]. In contrast,
the proposed spectral control approach, with approximate efficiency function, provides a solution for which the mean electric power, $\bar{P}_{e,\kappa|\mu}^{\dagger}$, lies within a few percents of the true optimal power value. Unlike other studies [11–14], the controller’s performance, optimised taking $\mu$ into account, is reasonably close to the “best-case scenario” $\mu\bar{P}_{id}^{*}$ (in fact, almost identical in sea states 1 and 4). Results also clearly evidence the interest of the proposed spectral control approach, with respect to an optimised linear damper, even in sea states 1 and 4, where the waves are close to the WEC resonant frequency.

Figure 3 shows the effect of a non-ideal PTO efficiency, upon the WEC optimal trajectory (in this case, represented by the WEC angular velocity $\dot{\theta}$, see Section 3.2), control torque, and instantaneous power, in Sea State 3. In the non-ideal case, $\mu$ is set to 0.7 as in the WECCCOMP.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Excitation torque, optimal angular velocity, control torque and instantaneous electric power, with $\mu = 1$ and $\mu = 0.7$.</figcaption></figure>
Although the two optimal velocities and control inputs look qualitatively similar<sup>6</sup>, the optimal electric power is fundamentally modified by the non-ideal PTO efficiency. As could be expected, very little reactive power must be employed with μ = 0.7; furthermore, the peaks in instantaneous electric power are reduced by a factor of at least 2 (while the average electric power still exceeds 60% of the ideal one).

## 4.2. The RM3 Device

Results for the RM3 device (see Section 3.3) are illustrated in a way similar to those of Section 4.1, in Figures 4 and 5. In Figure 4, the two bounds $\bar{P}_{e,\kappa}^*|_{\mu}$ and $\bar{P}_{e,\kappa}^\dagger|_{\mu}$ are close enough to indicate that the solution obtained is nearly optimal.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Normalised mean electric power for the RM3 device in four JONSWAP sea states with $\gamma = 3$, over a range of efficiency parameter values.</figcaption></figure>

The RM3 resonant period of the device is of the order of 4.5s, which is shorter than typical waves encountered at the device's design location<sup>7</sup>. As a consequence, from Figure 4, it can be seen that for shorter $T_p$, less reactive power is used to achieve optimal power absorption, thus making the simple linear damper perform reasonably well with respect to the reactive control strategy. In contrast, when waves are away from the WEC resonant frequencies, more reactive power must be employed, and the control strategy is significantly affected by the non-ideal efficiency; see, in particular, $\bar{P}_{e,mism.}$. However, the efficiency-aware spectral control strongly mitigates such adverse effects, bringing $\bar{P}_{e,\kappa}^\dagger$ close to the hypothetical best-case scenario $\mu \bar{P}_{id}^*$.

Finally, similarly to Figure 3, the effect of a non-ideal efficiency upon the control strategy is examined more closely in Figure 5, for the JONSWAP sea state with $T_p = 8$ s. While the ideal PTO

<footnote n="6"><sup>6</sup> In particular, in both cases the optimal velocity is in phase with the excitation torque, which is a classical WEC optimal control result [1]</footnote>

<footnote n="7"><sup>7</sup> This is typically the case for heaving point absorbers: their dimensions give them a relatively short resonant period, but reactive control can artificially make them resonate in longer wave periods.</footnote>
yields large reactive power flows, taking into account a non-ideal efficiency $\mu = 0.7$ eliminates most of the reactive power, and spectacularly reduces the peak power amplitude. With the non-ideal PTO efficiency, the optimal relative heave velocity $\dot{x}$ (see Section 3.3) and control force $u$ exhibit visible higher-frequency oscillations. Finally, note that the relative position constraint is successfully implemented within the spectral control strategy, which would not be possible with a simpler control approach, such as a PI controller.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Excitation force, optimal relative velocity, relative position, control force, and instantaneous electric power, with $\mu = 1$ and $\mu = 0.7$, $T_p = 8$ s.</figcaption></figure>

**5. Conclusions**

This work shows how the spectral control approach in [5] can be extended to handle more general objective functions, of the form outlined in Equation (2). The efficiency of the proposed technique allows optimal control results to be calculated at once in long wave signals, comprising more than
20 wave periods. Such calculations permit a theoretical analysis of WEC optimal control results, under a variety of modelling and control choices. Furthermore, as in [9], they can provide a robust point of comparison to assess the performance of more practical control approaches, typically implemented in a receding-horizon configuration. In fact, the authors have used the optimal spectral control approach proposed in this work, to assess the receding-horizon MPC solution, proposed by their institution for the WECCCOMP. This comparison will be the subject of a future paper.

A particularly important application of the proposed spectral method consists of WEC control with non-ideal PTO efficiency. In this study, a discontinuous efficiency function, given in Equation (16), is considered. A smooth, hyperbolic tangent approximation for the efficiency function allows the calculation of a near-optimal solution, which also provides upper and lower bounds to the exact optimal electric power.

In view of the results of Figures 2 and 4, concerning two different WECs representative of heaving point absorbers (both devices have resonant periods shorter than their design waves), a non-ideal PTO efficiency may not be as detrimental as suggested in previous studies [11,12,14,15], if the control strategy is able to take the non-ideal efficiency into account, which is the case for the proposed optimal spectral approach. In fact, with realistic efficiency values in the range of $0.7 \le \mu \le 0.9$ [11], it does not seem unreasonable to approximate the optimal electric power achievable with $\mu \bar{P}_{id}^*$, minus 5 to 15% (reminding that $\mu \bar{P}_{id}^*$ is merely a theoretical limit, that no control strategy can exceed if the PTO has an efficiency factor of $\mu$).

The fact that a non-ideal PTO efficiency necessitates a control strategy with less reactive power, already evidenced in previous studies [11,12], is clearly corroborated by the results of Figures 3 and 5. This is to be expected, irrespective to the WEC technology and the sea states considered. However, the extent of power production degradation when a control strategy neglects the non-ideal PTO characteristics will be, in general, system- and wave climate-dependent. Note, also, that the results of this study depend, to some extent, on the specific model adopted for the PTO efficiency. Other efficiency functions could be investigated, such as those in [10,13], which model the fact that the PTO efficiency drops for small load factors.

The computational speed of the proposed spectral control technique is not a central aspect of the present work; however, it was found that optimal solutions, for irregular wave signals with more than 5 wave periods, could be obtained in less than 1 second, with a Matlab implementation running on a 2.60 GHz, 4-core Intel® processor. A real-time, receding-horizon implementation could be considered in future work, either by employing windowing functions to make the receding-horizon control problem periodic, in the spirit of the works in [8,9], or through the use of non-periodic spectral basis functions [6].

**Author Contributions:** Conceptualising, A.M.; software, A.M.; writing—original draft preparation, A.M.; writing—review and editing, A.M. and P.T.; supervision, P.T.; project administration, P.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

**References**

1. Ringwood, J.V.; Bacelli, G.; Fusco, F. Energy-maximizing control of wave-energy converters: The development of control system technology to optimize their operation. *IEEE Control Syst. Mag.* **2014**, *34*, 30–55.

2. Faedo, N.; Olaya, S.; Ringwood, J.V. Optimal control, MPC and MPC-like algorithms for wave energy systems: An overview. *IFAC J. Syst. Control* **2017**, *1*, 37–56. [CrossRef]

3. Bacelli, G.; Ringwood, J.V.; Gilloteaux, J.C. A control system for a self-reacting point absorber wave energy converter subject to constraints. *IFAC Proc. Vol.* **2011**, *44*, 11387–11392. [CrossRef]

4. Bacelli, G.; Ringwood, J.V. Nonlinear optimal wave energy converter control with application to a flap-type device. *IFAC Proc. Vol.* **2014**, *47*, 7696–7701. [CrossRef]
5. Mérigaud, A.; Ringwood, J.V. Improving the computational performance of nonlinear pseudospectral control of wave energy converters. *IEEE Trans. Sustain. Energy* **2017**, *9*, 1419–1426. [CrossRef]

6. Genest, R.; Ringwood, J.V. Receding horizon pseudospectral control for energy maximization with application to wave energy devices. *IEEE Trans. Control Syst. Technol.* **2016**, *25*, 29–38. [CrossRef]

7. Boyd, J.P. *Chebyshev and Fourier Spectral Methods*; Courier Corporation: North Chelmsford, MA, USA, 2001.

8. Auger, C.; Merigaud, A.P.L.; Ringwood, J.V. Receding-horizon pseudo-spectral control of wave energy converters using periodic basis functions. *IEEE Trans. Sustain. Energy* **2018**, *10*, 1644–1652. [CrossRef]

9. Mérigaud, A.; Ringwood, J.V. Towards realistic nonlinear receding-horizon spectral control of wave energy converters. *Control Eng. Pract.* **2018**, *81*, 145–161.

10. Tedeschi, E.; Carraro, M.; Molinas, M.; Mattavelli, P. Effect of control strategies and power take-off efficiency on the power capture from sea waves. *IEEE Trans. Energy Convers.* **2011**, *26*, 1088–1098. [CrossRef]

11. Genest, R.; Bonnefoy, F.; Clément, A.H.; Babarit, A. Effect of non-ideal power take-off on the energy absorption of a reactively controlled one degree of freedom wave energy converter. *Appl. Ocean. Res.* **2014**, *48*, 236–243. [CrossRef]

12. Falcão, A.F.; Henriques, J.C. Effect of non-ideal power take-off efficiency on performance of single-and two-body reactively controlled wave energy converters. *J. Ocean. Eng. Mar. Energy* **2015**, *1*, 273–286.

13. Bacelli, G.; Genest, R.; Ringwood, J.V. Nonlinear control of flap-type wave energy converter with a non-ideal power take-off system. *Annu. Rev. Control* **2015**, *40*, 116–126. [CrossRef]

14. Sánchez, E.V.; Hansen, R.H.; Kramer, M.M. Control performance assessment and design of optimal control to harvest ocean energy. *IEEE J. Ocean. Eng.* **2014**, *40*, 15–26. [CrossRef]

15. Tom, N.M.; Madhi, F.; Yeung, R.W. Power-to-load balancing for heaving asymmetric wave-energy converters with nonideal power take-off. *Renew. Energy* **2019**, *131*, 1208–1225. [CrossRef]

16. Tona, P.; Nguyen, H.N.; Sabiron, G.; Creff, Y. An Efficiency-Aware Model Predictive Control Strategy for a Heaving Buoy Wave Energy Converter. In Proceedings of the European Wave and Tidal Energy Conference (EWTEC) 2015, Nantes, France, 6–11 September 2015.

17. Saupe, F.; Creff, Y.; Gilloteaux, J.C.; Bozonnet, P.; Tona, P. *Latching Control Strategies for a Heaving Buoy Wave Energy Generator in a Random Sea*; IFAC World Congress: Prague, Czech Republic, 2014; Volume 19, pp. 7710–7716.

18. Ringwood, J.V.; Ferri, F.; Ruehl, K.; Yu, Y.H.; Coe, R.G.; Bacelli, G.; Weber, J.; Kramer, M.M. A competition for WEC control systems. In Proceedings of the 12th European Wave and Tidal Energy ConferenceEuropean Wave and Tidal Energy Conference, Cork, Ireland, 27 August–1 September 2017.

19. Neary, V.S.; Lawson, M.; Previsic, M.; Copping, A.; Hallett, K.C.; LaBonte, A.; Rieks, J.; Murray, D. *Methodology for Design and Economic Analysis of Marine Energy Conversion (MEC) Technologies*; Technical Report; Sandia National Lab. (SNL-NM): Albuquerque, NM, USA, 2014.

20. Mérigaud, A. A Harmonic Balance Framework for the Numerical Simulation of Non-Linear Wave Energy Converter Models in Random Seas. Ph.D. Thesis, Maynooth University, Maynooth, Ireland, 2018.

21. Cummins, W. The impulse response function and ship motions. *Schiffstechnik* **1962**, *9*, 101–109.

22. Ogilvie, T. Recent progress toward the understanding and prediction of ship motions. In Proceedings of the Sixth Symposium on Naval Hydrodynamics, Bergen, Norway, 10–12 September 1964.

23. Hasselmann, K.; Barnett, T.; Bouws, E.; Carlson, H.; Cartwright, D.; Enke, K.; Ewing, J.; Gienapp, H.; Hasselmann, D.; Kruseman, P. *Measurements of Wind-Wave Growth and Swell Decay During the Joint North Sea Wave Project (JONSWAP)*; Technical Report; Deutsches Hydrographisches Institut: Hamburg, Germany, 1973.

24. Olaya, S. Contribution à la Modélisation Multi-Physique et au Contrôle Optimal d'un Générateur Houlomoteur: Application à un Système 'deux corps'. Ph.D. Thesis, Université de Bretagne occidentale, Brest, France, 2016. (In French)

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).