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

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.


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, 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ẋ, as highlighted by a number of authors [10][11][12][13][14][15][16][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 ≤ P e ≤ 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 ≤ 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 ≤ P a ≤ 0. PTO losses should be taken into account when designing the WEC control strategy [10][11][12][13][14][15][16][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 where α is some scalar function of x,ẋ and u (while the control calculations presented in [5] can only deal with cost functions in the form −uẋ, quadratic in the problem variables). P e , being modelled as a function of P a = −uẋ, 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.

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]): where the dependence in the time variable t has been omitted to alleviate notations, and • L represents a linear integro-differential operator, containing some of, or all the terms of the well-known Cummins' equation [21]: where I represents the system inertia; I ∞ 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 L is represented by a frequency-response mapping z(ω) = S − ω 2 (I + A r (ω)) − jω(C + B r (ω)), where the radiation frequency-dependent added mass A r (ω) and damping B r (ω) are related to I ∞ and k r through Ogilvie's relation [22].

Spectral Control Problem Formulation
Consider an optimal control problem of the following form, 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 α, 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 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 2n cos(ω n t) +ê 2n+1 sin(ω n t) where the frequencies ω n are defined harmonically with ∆ω = 2π T and ∀n ∈ {1...N}, ω n = n∆ω. 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 2n cos(ω n t) +x 2n+1 sin(ω n t) 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ê n , except forê 2 andê 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 ≥ 2N + 1, so that the maximisation problem becomes where, for the sake of conciseness, α(t m ) denotes the value of α when its arguments are evaluated at time t m .
Furthermore, the linear terms of L evaluated at times t 1 .
Finally, the objective function can be expressed as a function ofx, as follows, where 1 M×1 is the unit vector of size M × 1 and, with a slight abuse of notation, the functions α and n are applied component-wise to the vectors used as arguments, i.e., for example, Note that Equation (13) differs from the objective function in [5], as α is, in general, not a quadratic function of its arguments. The calculations which follow constitute a useful generalisation of the results in [5].

Explicit Calculation of the Objective Function Derivatives
In this section, it is assumed that α and n are continuous in all their arguments, and admit second-order derivatives. Under such an assumption, the first-and second-order derivatives off , with respect to the components ofx, 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 × M diagonal matrix with diagonal entries y(t 1 ), ...y(t M ). Then, the vector of 1st-order derivatives is given as with The matrix of second-order derivatives is given as with: In practice, at each iteration of the optimisation algorithm, the objective functionα, the nonlinear terms n, and their 1st-and 2nd-order derivatives, are evaluated along the trajectory given by x = Φx.
f is evaluated; v 1 ...v 3 are built in order to calculate ∂f ∂x as in Equation (14), and w 1 ...w 6 are built in order to calculate ∂ 2f ∂x 2 as in Equation (15).

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.
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ẋ − 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ẋ, 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 µ is, the less reactive control should be performed (i.e., |k| should be smaller).
Overall, all studies [10][11][12][13][14][15]17,18] are relatively pessimistic, regarding the mean electrical power P e which can be achieved with a non-ideal PTO system: a 10% decrease in efficiency results in a drop in P e of more than 50%, even when control parameters are tuned taking the PTO efficiency into account 1 .
Define h(P a ) the efficiency function, so that h(P a ) = P e P a . Consistently with Equation (16), h is equal to µ for positive P a , and 1/µ 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, where A = 1 2 (µ − 1/µ), B = 1 2 (µ + 1µ) and κ > 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. If the exact efficiency is considered, the function α in the optimisation problem (5) is α(ẋ, u) := −uẋ h(−uẋ). However, in practice, the problem is solved using the approximate function α κ (ẋ, u) := −uẋ h κ (−uẋ). For all (ẋ, u), α(ẋ, u) ≤ α κ (ẋ, u). Therefore, the optimal mean electric power, P * e,κ , obtained with α κ , is an upper bound to the optimal exact electric power, P * e , which would be obtained if Problem (1) could be solved using the exact efficiency. Furthermore, for all κ 1 ≤ κ 2 , for all (ẋ, u), α κ 2 (ẋ, u) ≤ α κ 1 (ẋ, u), which implies that P * e,κ decreases monotonically with κ, and has P * e as a limit when κ → ∞.
P * e can also be bounded from below: consider the solution x * κ , u * κ , obtained with the approximate objective function. x * κ , u * κ is necessarily sub-optimal, with respect to the exact objective function. Therefore, the exact mean electric power, effectively achieved with (x * κ , u * κ ), constitutes a lower bound to P * e , and is calculated as follows.
It should be noted, however, that the control strategies adopted, and the PTO efficiency models, differ across studies, so that accurate comparisons are difficult. Therefore the following inequality holds.
Even if the true optimal electrical power P * e cannot be calculated for the exact efficiency, by setting κ to a large enough value, it can be ensured that the electric power P † e,κ , effectively achieved with the solution found, lies within a prescribed percentage of P * e . Finally, let (x * , u * )| µ be the optimal solution for efficiency µ < 1, and P * e | µ (resp. P a | µ ) the mean electrical (resp. absorbed) power with the solution (x * , u * )| µ . It is easy to show 2 that P * e | µ < µP a | µ . Define P * id as the optimal power for an ideal PTO (absorbed and electric power are then identical). The solution (x * , u * )| µ is sub-optimal with respect to an ideal PTO, thus P a | µ < P * id . Therefore, the following simple inequality holds.
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 µ < 1, and the inequality becomes an equality: P * e | µ = µ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.

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 = 40s, at least equal to 20 typical wave periods for the sea states considered.
In each random realisation, the optimal solutionx * κ | µ , and power P * e,κ | µ , with µ = 1, 0.9, . . . , 0.5, are calculated using the approximate efficiency function h κ . Then, the exact power P † e,κ | µ , obtained witĥ x * κ | µ , is also calculated. For each pair of sea states (1-4, 2-5, 3-6), the appropriate range of values for κ differs (when typical values for absorbed power are larger, a smaller κ yields a satisfactory approximation). For each pair of sea states, κ is gradually increased, until bounds given by Equation (19) are sufficiently close. Suitable values for κ are indicated in Table 1. To see this, decompose P a into its passive and reactive components. 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 µ between 1 and 0.5. In each sea state, power values are normalised by P * id , the optimal power obtained with µ = 1. 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 P * id = ∞ ω=0 −ω|ê(ω)| 2 2 {z(ω)} dω. As per Equation (19), P † e,κ | µ and P * e,κ | µ together provide bounds to P * e | µ , the optimal power for the exact efficiency function. Also shown in Figure 2, the dotted lines labelled "µP * id " represent the upper bound of Equation (20), i.e., a hypothetical scenario where P * id would "only" be affected by a factor of µ. P mism. is the actual electric power which would be obtained, if the solutionx * id to the ideal control problem were followed by a WEC with imperfect PTO efficiency µ ≤ 1; in other words, P mism. shows the effect of a mismatch between the nominal efficiency (µ = 1) and the actual efficiency (1 > µ ≥ 0.5). Finally, P Bpto represents the power obtained using a simple linear PTO damping, with coefficient B pto optimised for each sea state.

The WECCCOMP Model
The importance of taking into account the non-ideal PTO efficiency in control calculations is evidenced by 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, P † e,κ | µ , lies within a few percents of the true optimal power value. Unlike other studies [11][12][13][14], the controller's performance, optimised taking µ into account, is reasonably close to the "best-case scenario" µ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θ, see Section 3.2), control torque, and instantaneous power, in Sea State 3. In the non-ideal case, µ is set to 0.7 as in the WECCCOMP. Although the two optimal velocities and control inputs look qualitatively similar 6 , 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).

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 P * e,κ | µ and P † e,κ | µ are close enough to indicate that the solution obtained is nearly optimal. 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 7 . 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, P e,mism. . However, the efficiency-aware spectral control strongly mitigates such adverse effects, bringing P † e,κ close to the hypothetical best-case scenario µ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 6 In particular, in both cases the optimal velocity is in phase with the excitation torque, which is a classical WEC optimal control result [1] 7 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. yields large reactive power flows, taking into account a non-ideal efficiency µ = 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ẋ (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 5. Excitation force, optimal relative velocity, relative position, control force, and instantaneous electric power, with µ = 1 and µ = 0.7, T p = 8 s.

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 ≤ µ ≤ 0.9 [11], it does not seem unreasonable to approximate the optimal electric power achievable with µP * id , minus 5 to 15% (reminding that µP * id is merely a theoretical limit, that no control strategy can exceed if the PTO has an efficiency factor of µ).
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].