## Abstract

This paper presents a solution to the optimal control problem of a three degrees-of-freedom (3DOF) wave energy converter (WEC). The three modes are the heave, pitch, and surge. The dynamic model is characterized by a coupling between the pitch and surge modes, while the heave is decoupled. The heave, however, excites the pitch motion through nonlinear parametric excitation in the pitch mode. This paper uses Fourier series (FS) as basis functions to approximate the states and the control. A simplified model is first used where the parametric excitation term is neglected and a closed-form solution for the optimal control is developed. For the parametrically excited case, a sequential quadratic programming approach is implemented to solve for the optimal control numerically. Numerical results show that the harvested energy from three modes is greater than three times the harvested energy from the heave mode alone. Moreover, the harvested energy using a control that accounts for the parametric excitation is significantly higher than the energy harvested when neglecting this nonlinear parametric excitation term.

## Introduction

Despite its potential, wave energy extraction technology is not yet mature due to few challenges. The Department of Energy reported in 2012 that the world wave installed capacity is less than 5 MW [1], mostly produced on a precommercial basis. One of the main challenges is the design of an effective control strategy that enables the harvesting of higher levels of energy. Another challenge is the limitation of the existing theory and practice. Energy harvesting is about the wave energy converter (WEC) motion and control. The fundamental theory that is being used to describe this motion is inherited from ship hydrodynamics. Ships are usually designed so as to not move in the vertical (heave) direction or rotate in the pitch or roll directions. Hence most theories, designs, and control algorithms are based on the assumption of small motion in these directions, which yield the time-invariant linear dynamic models. WEC devices, however, need to heave and/or pitch in order to harvest more energy. This paper uses a dynamic model that accounts for a second-order kinematic nonlinearity and shows that it can be handled as a parametric excited time-varying linear system. A pseudo-spectral (PS) approach is then used to search for the optimal control.

Usually, a linear system is used to model the dynamics of WECs, as detailed in several references in the literature, e.g., see Refs. [2–5]. Using a linear WEC model, significant developments have been made on the control design of single degree-of-freedom WECs. The latching control [6] is one classical method for WEC control. A model predictive control is developed in multiple references, e.g., see Ref. [7], that takes hardware limitation into account. Recently, a multi resonant control was introduced in Ref. [8] that is a time domain implementation of the complex conjugate control [3].

Several references have motivated the use of a multiple degrees-of-freedom (multi-DOF) WEC as opposed to a single-mode WEC. Evans [9] extended the results of two-dimensional WECs to bodies in channels and accounts for the body orientation on the energy harvesting. In fact, Ref. [10] points out that the power that can be extracted from a mode that is antisymmetric to the wave (such as pitch and surge) is twice as much as that can be extracted from a mode that is symmetric (such as heave). One of the references that recently studied the pitch–surge power conversion is Ref. [11]. Yavuz [11] models the pitch–surge motions assuming no heave motion; hence, there is no effect from the heave motion on the pitch–surge power conversion. The mathematical model used in Ref. [11] for the motions in these 2DOF WECs is coupled through mass and damping only; there is no coupling in the stiffness.

However, it has been observed that floating structures can be subject to parametric instability arising from variations of the pitch restoring coefficient [12]. Villegas and van der Schaaf [12] describe an experimental heaving buoy for which the parametric excitation causes the pitch motion to grow resulting in instability. A harmonic balance approach is implemented to cancel this parametric resonance and results of tank experiments are presented [12]. As will be detailed in this paper, the equations of motion for a 3DOF WEC have a second-order term that causes the heave motion to parametric excite the pitch mode; and the pitch and surge motions are coupled. For relatively large heave motions, which would be needed for higher energy harvesting, it is not possible to neglect this parametric excitation term. Rather, the controller should be designed to leverage this nonlinear phenomenon for optimum energy harvesting.

In designing the optimal control, it is important to note that the WEC buoy is excited by waves which are stochastic in nature. The dynamic programming approach is one method that can be implemented to search for the optimal control in such a stochastic system [13]. The dynamic programming approach, however, is computationally expensive and suffers the curse of dimensionality. Indirect methods for optimal control seek a solution to the closed system of necessary conditions of optimality [14]. Direct methods for optimal control discretize the control problem and then apply nonlinear programming (NLP) techniques to the resulting finite-dimensional optimization problem [15]. Due to control parameterization, direct methods provide only approximate solutions.

This paper presents a control method for nonlinear 3DOF WECs. The contribution of this paper is: (1) a 3DOF model that accounts for the nonlinear coupling between the heave and pitch modes is used and (2) a PS control [16] on this 3DOF WEC model is implemented. In this paper, a Fourier series (FS) basis is used to approximate the control and the states. The objective function maximizes the extracted energy. Two cases will be investigated; the first case assumes a linear dynamic model where no parametric excitation is modeled. The second case uses a model that has parametric excitation. The obtained results are compared to the results obtained using a nonlinear model predictive control (NMPC) approach for comparison.

## Mathematical Model

The dynamic model of a 3DOF WEC is here briefed. Figure 1 shows a typical geometry for a cylindrical buoy, where the surge displacement is *d*_{1}, the heave displacement is *d*_{3}, and the pitch rotation angle is *θ*_{5}. The cylindrical buoy has a radius *R* and a mass *m*. The height between the center of gravity *G* and the bottom of the buoy is *h*. *B* denotes the center of buoyancy.

*i*th mode. The radiation forces can be expressed as

where * denotes the convolution operator, and *h _{ij}* are the radiation impulses.

*M*

_{total}] is summation of the mass of the rigid body [

*M*], and the added mass at infinite frequency [

*M*]

_{∞}*C*] is also expressed as the summation of viscous damping and radiation damping

*K*] is

*K*

_{11}=

*K*

_{moor}is the mooring stiffness in surge direction, and

The time-varying part of the stiffness, *K _{p}*(

*t*), can be considered as a parametrically excited term [19]. As can be seen from the equations of motion above, the heave motion excites the pitch motion; also pitch and surge modes are coupled.

## Pseudospectral Optimal Control

### System Approximation Using Fourier Series.

*ϕ*(

_{k}*t*). For the WEC problem, and due to the periodicity nature of the wave, it is intuitive to select a Fourier series to be the basis functions. A truncated Fourier series that has zero mean is used with

*N*terms. Since we have both sine and cosine functions,

*N*is an even number and is equal to twice the number of cosine (or sine) functions in the Fourier series. The vector of basis functions is

*w*

_{0}= 2

*π*/

*T*

_{end}is the fundamental frequency and

*T*

_{end}represents the total simulation time. The states and the controller can be approximated using Φ(

*t*) as follows:

*i*is the state index and

*j*is the control index. In Eq. (9), $xikc/xiks$ denotes the

*k*th coefficient of cosine/sine term of basis function for

*i*th state. In Eq. (10), $ujkc/ujks$ denotes the

*k*th coefficient of cosine/sine term of basis function for

*j*th control. The Fourier coefficients (or weight vectors) are grouped as follows:

*v*(

_{s}*t*), the velocity of pitch rotation

*v*(

_{p}*t*), the controller in surge direction

*u*(

_{s}*t*), and the controller in pitch direction

*u*(

_{p}*t*). So, we have

*i*= 1, 2

*j*= 1, 2

*D*∈

_{ϕ}*R*

^{N}^{×}

*, the matrix is block diagonal, where each block $D\varphi k$ can be expressed as*

^{N}*D*is invertible, and its inverse is the matrix used to compute the integration of a state. The integration matrix is still block diagonal. Each block of integration matrix can be written as

_{ϕ}*f*can be expressed as

_{h}*m*denotes the added mass at infinite frequency. The matrix

_{∞}*G*∈

*R*

^{N}^{×}

*is block diagonal, the*

^{N}*k*th block is

*r*

_{1}and

*r*

_{2}are residuals due to the approximation. In the Galerkin method, the residuals are orthogonal to the basis function. That is

where *j* = 1,…, *N*. This implies that the residuals vanish at all the collocation points.

### Control Optimization for Linear Time-Invariant Wave Energy Converters.

*J*defined as

*F*are given as a block diagonal matrix

_{ij}*i*,

*j*= 1, 2

where *A _{ij}*(

*kω*

_{0}) denotes the added mass at

*k*th frequency of

*i*th mode due to

*j*th mode.

### Control Optimization for Linear Time-Varying Wave Energy Converters.

*K*

_{55}has a time-varying term

*K*(

_{p}*t*). This is usually referred to as parametric excitation since the heave motion excites the pitch motion, which is coupled with the surge. Then the system becomes a time-varying system. This affects the second constraint Eq. (29), which becomes as follows in a linear time-varying WEC:

*r*

_{2}residual can be expressed as

*r*

_{1}residual remains the same. The residuals can be discretized at collocation points

*= Φ(*

_{j}*t*), the nodes

_{j}*t*are uniformly spaced between 0 and

_{j}*T*

_{end}

The objective function will be the same as the time-invariant system in Eq. (33). To solve for the control, an NLP approach will be implemented. The NLP is then to minimize Eq. (33), subject to Eqs. (40) and (41), which have a total of 2*N* equality constraints. Figure 2 shows a framework for the control logic and the simulation environment is used to test the proposed control. The switch is on when a new control update is needed.

## Numerical Results

In all the numerical test cases, a buoy is assumed of cylindrical shape with a tapered end as depicted in Fig. 3, which is similar to the buoy of the Sandia National Lab [22]. The radius of the buoy is 0.8604 m in the cylindrical part, and its cross section area is 2.3257 m^{2}. The mass of the buoy is 858.3987 kg.

Both regular and irregular waves are tested. The regular wave has a wave amplitude of 0.4 m. The irregular wave has a Bretschneider spectrum, with a significant height of 0.5 m and a peak period of 1.5708 s. The Bretschneider spectrum is realized using 1260 frequencies. The step for updating the control is set to 0.01 s. The heave motion of the buoy is constrained to be within ±0.2 m; this way, the water line will not reach the conical section of the buoy and the dynamic model for cylindrical buoys, adopted in this paper, remains valid.

The comparison is made between the system with and without parametric excitation. The details for the velocity of surge and pitch motions are presented. The velocities are plotted along with the excitation force to analyze the system response.

The system parameters and the control parameters are listed in the table below.

### No Parametric Excitation.

This section presents numerical results for the linear time-invariant WEC. A regular wave is used in this simulation. The number of basis functions is *N* = 80 in this case. The energy captured by surge, heave, and pitch is shown in Fig. 4, respectively. The total absorbed energy is around 3.5 times the energy absorbed from only heave mode. It is to be noted that this amplification factor is obtained when we take into account the viscous damping. If viscous damping is removed, then this amplification factor reduces to only three. The control force is presented in Fig. 5 while the velocity of surge and pitch modes, respectively, is shown in Figs. 6 and 7. In these two figures, the excitation force is scaled to show the resonance between the velocity and the excitation force, which indicates that optimality is achieved. In these figures, both the true velocities and the approximated velocities are plotted.

System parameters | |

Mass (kg) | 858.3987 |

Moment of inertia (kg/m^{2}) | 84 |

Hydrostatic stiffness (N/m) | 2.3981 × 10^{4} |

Mooring stiffness (N/m) | 126.6462 |

Surge viscous damping coefficient (N s/m) | 156.4364 |

Heave viscous damping coefficient (N s/m) | 253.9528 |

Pitch viscous damping coefficient (N.m.s/r) | 37.0737 |

Control parameters | |

Number of collocation points | TBD |

Maximum surge displacement (m) | |

(constrained) | 1 |

Maximum pitch angle (r) | |

(constrained) | 1 |

System parameters | |

Mass (kg) | 858.3987 |

Moment of inertia (kg/m^{2}) | 84 |

Hydrostatic stiffness (N/m) | 2.3981 × 10^{4} |

Mooring stiffness (N/m) | 126.6462 |

Surge viscous damping coefficient (N s/m) | 156.4364 |

Heave viscous damping coefficient (N s/m) | 253.9528 |

Pitch viscous damping coefficient (N.m.s/r) | 37.0737 |

Control parameters | |

Number of collocation points | TBD |

Maximum surge displacement (m) | |

(constrained) | 1 |

Maximum pitch angle (r) | |

(constrained) | 1 |

### Parametric Excited Surge–Pitch Motions.

This section presents numerical results for the linear time-varying WEC. A regular wave is used in this simulation. A nonlinear programming approach is implemented for computing the controller numerically. The number of basis functions is *N* = 80. The energy captured by the three modes is shown in Fig. 8. The energy extracted from surge mode is 3.661 × 10^{4} J, from the heave mode is 1.401 × 10^{4} J, and from the pitch mode is 1.701 × 10^{4} J within 30 s. The energy captured in the three modes is about 4.8 times the energy captured by heave mode only. Compared to the nonparametric excited system presented in Sec. 4.1, this parametrically excited system captures much more energy, about 1.4 times that of the nonparametric excited system. The controller is shown in Fig. 9, and the velocity for both pitch and surge modes are shown in Figs. 10 and 11, respectively.

### Parametric Excited Wave Energy Converter in Bretschneider Wave.

A more realistic case is simulated using an irregular wave. A Bretschneider spectrum of 0.5 m significant height, 1.5708 s peak period, and random initial phase for different frequencies is used in this section. The wave amplitude versus frequency is shown in Fig. 12. The number of basis functions is *N* = 120 in this case. Table 1 shows the system and control parameters used in this case.

Figure 13 shows the surge, heave, and pitch energy captured when we use a Bretschneider wave. The pitch energy is the highest energy compared to the other motions. There is some transient response initially, but the energy keeps growing at steady-state. The total energy is around 6.9 times the energy from heave mode only. The surge position and pitch rotation are shown in Figs. 14 and 15, respectively. The velocity of surge motion and the angular velocity of pitch rotation are shown in Figs. 16 and 17. Because of the parametric excitation, the optimal velocity of the pitch motion is not in resonance with the excitation force. The control forces along the pitch and surge modes are plotted in Fig. 18. It is noted here that this paper presents a method for computing the control; it is beyond the scope of this paper how the power take off unit can generate this control force.

### Constraint Optimization.

As can be seen from the results in Sec. 4.3, the pitch, heave, and/or surge motions might get larger than what is feasible in terms of the device hardware and operational limitations. It is straightforward to accommodate constraints on the states when implementing this pseudo-spectral approach for optimizing the control. Consider the NLP described at the end of Sec. 3.3. The state constraints can be appended to the constraints in Eqs. (40) and (41), and the problem is solved using the same approach as detailed above, but with more constraints. In this section, results are presented for the same buoy above, but limiting the surge displacement to 1 m, and the pitch angle to 1 rad. Figure 19 shows the extracted energy in each of the 3DOF for a parametric exited WEC in a Bretschneider wave. The energy harvested from a 3DOF WEC, in this constrained-motion case, is 4.6 times that harvested from a heave-only WEC. Figure 20 shows the pitch motion as a demonstration for the effectiveness of the controller in satisfying the motion constrains.

### Comparison Between Pseudo-Spectral and Nonlinear Model Predictive Control Methods.

This section compares the proposed FS pseudo-spectral method to an NMPC [18] approach in terms of optimizing the objective function. The surge and pitch motions are constrained to be smaller than 1 m and 1 rad, respectively. First, a regular wave, which has a wave elevation of 0.4 m, is used to test both controllers for a nonparametric excited WEC, and the harvested energy in each case is shown in Fig. 21. Figure 22, on the other hand, shows the harvested energy for both controllers for a parametric excited WEC in a regular wave. Using an irregular wave with a Bretschneider spectrum of a significant wave height of 1.2 m, both controllers were tested on a parametric excited WEC, and the harvested energy is plotted in Fig. 23. In all the cases above, the pseudo-spectral method was able to harvest energy at about the same level or higher than the NMPC.

## Discussion

The pseudo-spectral approach implemented in this paper is essentially an approximate method; hence, there is an error. This error materializes as an error in computing the optimal control. In all the previous figures, the results are accurate since they are a simulation for the computed control with the dynamic model. Yet the performance might not be optimal due to the approximations in the pseudo-spectral approach. The approximation error varies depending on the problem. For instance, the approximation error in the constraint problem is different from that in the unconstraint problem. To see that, consider the results presented in Fig. 21. The energy captured from surge and pitch in the constrained motion case is 3.724 × 10^{4} J and it is higher than the energy captured in the unconstrained motion case (3.455 × 10^{4} J, see Fig. 8), which is not expected. The reason for that is the approximation error in pseudo-spectral method. In this case, the approximation error in the unconstrained motion case is higher and hence the solution is not as close to the optimal solution as that of the constrained motion case. For verification, the error in each state between the simulated true values and the approximated values is computed. Figure 24 shows this error for one state (the surge velocity), for different surge constrain values (Constraint in Fig. 24 is the maximum allowed surge). As can be seen from Fig. 24, the unconstrained case has higher error than the constrained ones, in this case. It is noted here that this conclusion would vary from one test case to another. In this test case for instance, the surge reaches a maximum of about 0.5 m when it is unconstraint. In other test cases for other sea states, the conclusion might vary.

## Conclusion

This paper presented an implementation for the pseudo-spectral method to search for the optimal control of a three degrees-of-freedom wave energy converter, characterized by a time-varying linear dynamic model. The numerical results show that using a dynamic model that accounts for the parametric excitation can be exploited, in designing the control, to increase energy harvesting. The results also confirm that using a three-degrees-of-freedom power take off unit results in harvesting energy significantly more than three times the amount of energy that can be harvested with a one-degree-of-freedom heave power take off unit. This pseudo-spectral method is computationally efficient and can be used in real-time implementations. In future work, the authors plan to design nonlinear controls for nonlinear WEC models that account for nonlinear hydrodynamics.

## Acknowledgment

Sandia National Laboratories is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000. The resilience nonlinear control project at Sandia National Labs is supported by The Department of Energy (DOE), Office of Energy Efficiency and Renewable Energy (EERE), Wind and Water Power Technologies Office (WWPTO).

## Funding Data

Sandia National Laboratories (Contract No. 1514349).