Control of Wave Energy Converters with Discrete Displacement Hydraulic Power Take-Off Units

Abstract

The control of ocean Wave Energy Converters (WECs) impacts the harvested energy. Several control methods have been developed over the past few decades to maximize the harvested energy. Many of these methods were developed based on an unconstrained dynamic model assuming an ideal power take-off (PTO) unit. This study presents numerical tests and comparisons of a few recently developed control methods. The testing is conducted using a numerical simulator that simulates a hydraulic PTO. The PTO imposes constraints on the maximum attainable control force and maximum stroke. In addition, the PTO has its own dynamics, which may impact the performance of some control strategies.

1. Introduction and Background

The control of ocean WECs has received a great deal of attention over the past several decades. In recent years in particular there have been significant developments with respect to maximizing the harvested energy from WECs from a control system analysis and design perspective. Many of these control methods were proposed for ideal conditions in the absence of stroke or force limitations, and assuming ideal power take-off (PTO) units. This paper presents comparisons between some recently developed control methods; these simulations include a model for a hydraulic PTO, which consequently imposes constraints on the displacement of the buoy and on the maximum possible control force. These simulations also highlight some insight regarding the needed reactive power for some of the discussed control methods. Several earlier controllers were developed in [1,2,3] for WECs with hydraulic PTOs. A hydraulic system was validated using AMEsim in [4]. The control methods discussed in this paper are tested using hydraulic PTOs for the first time in this paper.

This section presents a review of hydraulic PTO units. Figure 1 is a general layout for a typical hydraulic PTO. The hydraulic system is composed of the actuator, the valve, the accumulators and the motor. The motion of the buoy will compress/decompress the chamber of the actuator and transfer the wave power to the hydraulic system. All the hydraulic systems can be categorized into three main groups: the constant-pressure, the variable-pressure, and the constant–variable pressure hydraulic systems [5,6].

1.1. Constant Pressure Configuration

The first configuration is constructed with a low-pressure (LP) accumulator and a high-pressure (HP) accumulator. This type of hydraulic system can be achieved with a simple mechanism, and the control level is low.

The typical configuration of a constant-pressure hydraulic system is presented in detail in [7,8], using phase control. Control of the constant-pressure hydraulic system is achieved by implementing auxiliary accumulators in [3]. The latching and declutching controls are demonstrated in [9] using a constant-pressure hydraulic system. Additionally, a declutching control is presented in [10] for controlling a hydraulic PTO by switching on and off using a by-pass valve. The method is also tested with the SEAREV WEC with an even higher energy absorption. A detailed image of a single acting hydraulic PTO system with the phase control is presented in [11,12]. The hydraulic system implemented in SEAREV is presented in [13]. In [14], a novel model of the hydraulic PTO of the Pelamis WEC is developed, with the ability to apply reactive power for impedance matching. In [15], a double-action WEC with an inverse pendulum is proposed. The authors of [15] report that a double-action PTO can supply the output power in each wave period without large instantaneous fluctuating power. A double-acting hydraulic cylinder array is developed in [16], where the model is found to be adaptive to different sea states to achieve higher energy extraction. The authors of [17] present the optimization of a hydraulic PTO of a WEC for an irregular wave where optimal damping is achieved by altering the displacement of the variable-displacement hydraulic motor. The authors of [18] present the design and testing of a hybrid WEC that obtains higher energy absorption than a single oscillating body with a hydraulic PTO. A discrete-displacement hydraulic PTO system is studied in [19] for the Wavestar WEC. An energy conversion efficiency of $70\%$ is achieved. Additionally, adjustment of the force applied by the PTO is accomplished through implementing multiple chambers.

1.2. Variable Pressure

A variable-pressure hydraulic system is suggested in [20,21,22]. In this situation, the piston is connected directly to a hydraulic motor. This system can achieve better controllability, but the fluctuation of the output power is not negligible. Two hydraulic PTO systems are compared in [23], where constant-pressure hydraulic PTO and variable-pressure hydraulic PTO systems are compared. It was shown that a variable-pressure hydraulic PTO system would have a higher efficiency. The variable-pressure approach was also investigated in [24], where the hydraulic motor is used in order to remove the accumulator and control the output using the generator directly. A comparison between a constant-pressure system and a variable-pressure system was conducted in [25]; validation was conducted using AMEsim and demonstrated a good agreement. Power smoothing was achieved in [26] by means of energy storage.

1.3. Variable–Constant Pressure

The variable–constant pressure hydraulic system is constructed with two parts: the variable pressure part and the constant pressure part. The variable pressure part is accomplished by a hydraulic transformer. A generic oil hydraulic PTO system, applied to different WECs, is introduced in [27]. The authors of [28] developed a PID controller, with the reactive power supplied by the hydraulic transformer (working as a pump). A suboptimal control is suggested in [28] for practical implementation in terms of the efficiency of the PTO.

2. WEC Dynamics

In this section, the WEC dynamic model used in this paper is briefed. In this paper, the floater used in the simulations is the Wavestar absorber [9]. The floater has a single degree of freedom motion which is the pitch rotation. The geometry of the proposed absorber is depicted in Figure 2.

The WEC dynamic model can be described based on linear wave theory by Equations (1) to (9):

$$J_{rigid}\ddot{\theta }={\tau }_{ex}+{\tau }_{res}+{\tau }_{rad}-{\tau }_G-{\tau }_{PTO}$$

(1)

where $J_{rigid}$ is the moment of inertia of the rigid body. $\theta$ is the pitch rotation of the floater. ${\tau }_{ex}$ is the wave excitation torque acting on the buoy, ${\tau }_{res}$ is the restoring momentum, ${\tau }_{rad}$ is the radiation torque, and ${\tau }_G$ is the torque caused by the gravity. The PTO torque is ${\tau }_{PTO}$ which is applied by the hydraulic cylinder. The equation of motion can be further expanded as:

$$\ddot{\theta }=\frac{1}{J_{rigid}+J_{\infty }}({\tau }_{ex}-{\tau }_{PTO}-K_{res}\theta -h_r\ast \dot{\theta })$$

(2)

where $J_{\infty }$ is the moment of the added mass at infinite frequency, $K_{res}$ is the coefficient of the hydro-static restoring torque, and $h_r$ is the radiation impulse response function. In Equation (2), the radiation torque is expanded as:

$$\begin{array}{cc}\hfill {\tau }_{rad}& =-J_{\infty }\ddot{\theta }-{\tilde{\tau }}_{rad}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\tilde{\tau }}_{rad}& =h_r\ast \dot{\theta }\hfill \end{array}$$

(4)

The ∗ operation is the convolution between the impulse response function and the angular velocity $\dot{\theta }$ which can be approximated by a state space model as:

$$\begin{array}{cc}\hfill {\dot{x}}_r& =A_r{x}_r+B_r\dot{\theta }\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\tilde{\tau }}_{rad}& =C_r{x}_r+D_r\dot{\theta }\hfill \end{array}$$

(6)

where $A_r$, $B_r$, $C_r$, and $D_r$ are the radiation matrices which are identified from the radiation impulse response function. The excitation torque can be expressed by the convolution between the impulse response function and the wave elevation:

$${\tau }_{ex}=h_{ex}\ast \eta$$

(7)

Hence, the convolution can also be approximated by a state space model as:

$$\begin{array}{cc}\hfill {\dot{x}}_e& =A_e{x}_e+B_e\eta \hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\tau }_{ex}& =C_e{x}_e\hfill \end{array}$$

(9)

where $A_e$, $B_e$, and $C_e$ are the excitation matrices which are identified from the excitation impulse response function. The parameters of the floater are listed in

Table 1. Model parameters for the Wavestar.

Symbol	Value	Unit
$J_{tot}$	$3.8\times {10}^6$	kg m${}^2$
$K_{res}$	$14\times {10}^6$	Nm/rad
The transfer function $H_r(s)$
$(b_0,b_1,\dots ,b_5)$	*
$(a_0,a_1,\dots ,a_5)$	**
The transfer function $H_{ex}(s)$
$(b_0,b_1)$	***
$(a_0,a_1,\dots ,a_4)$	****

* = $(0.0001,0.0144,0.624,8.16,13.1,1.44)\times {10}^6$; ** = $(0.001,0.0906,1.67,6.31,13.3,9.18)$; *** = $(5.4,270)\times {10}^4$; **** = $(0.036,0.39,1.5,2.6,1.6)$.

. The viscous damping is not considered in the proposed dynamic model because it is assumed to be negligible based on linear wave theory. In this designed study the extreme wave motion will not be achieved due to the limited capacity of the PTO unit. As a result, the small wave assumption can be held. The frequency response of the proposed WEC dynamic model without control is shown in Figure 3.

3. The Hydraulic PTO System

In this paper, the Discrete Displacement Cylinder (DDC) hydraulic system is used to apply the PTO torque. A simplified illustration for this system is shown in Figure 4. More details about the DDC hydraulic system can be found in [19].

As shown in Figure 4, the DDC hydraulic system is mainly composed of the actuator/cylinder, the manifold valves, the manifold accumulators, and the generator. The PTO torque is computed in the Equation (10) as the product of the cylinder force and the moment arm:

$${\tau }_{PTO}=F_c{l}_1$$

(10)

where the moment arm can be expressed by Equation (11):

$$\begin{array}{cc}\hfill l_1& =\frac{l_2{l}_{3}sin(\theta -{\alpha }_0)}{x_c+l_4}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill x_c& =-l_4+\sqrt{-2l_2{l}_{3}cos(\theta -{\alpha }_0)+(l_2^2+l_3^2)}\hfill \end{array}$$

(12)

3.1. The Hydraulic Cylinder

The actuator force $F_c$ is generated by the hydraulic cylinder and can be computed by Equations (13) and (14):

$$\begin{array}{cc}\hfill {\tilde{F}}_c& =-p_{A1}A1+p_{A2}A2-p_{A3}A3+p_{A4}A4\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill F_c& ={\tilde{F}}_c-F_{fric}\hfill \end{array}$$

(14)

where $p_{Ai}$ is the pressure of the ith chamber and $Ai$ is the area of the piston. $F_{fric}$ is the cylinder friction force. The dynamics of the chamber pressure can be described by flow continuity Equations (15) to (18):

$$\begin{array}{cc}\hfill {\dot{p}}_{A1}& =\frac{\beta (p_{A1})}{A1(x_{c,max}-x_c)+V_{0,A1}}(Q_{A1}-v_{c}A1)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\dot{p}}_{A2}& =\frac{\beta (p_{A2})}{A2x_c+V_{0,A2}}(Q_{A2}+v_{c}A2)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\dot{p}}_{A3}& =\frac{\beta (p_{A3})}{A3(x_{c,max}-x_c)+V_{0,A3}}(Q_{A3}-v_{c}A3)\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\dot{p}}_{A4}& =\frac{\beta (p_{A4})}{A4x_c+V_{0,A4}}(Q_{A4}+v_{c}A4)\hfill \end{array}$$

(18)

where $V_{0,A1}$, $V_{0,A2}$, $V_{0,A3}$, and $V_{0,A4}$ are the volumes of the connecting hoses of different chambers. $x_{c,max}$ is the maximum stroke of the cylinder. $x_c$ and $v_c$ are the position and velocity of the piston, respectively, which are defined positive down. $\beta (p_{Ai})$ is the effective bulk modulus of the fluid based on different pressures, and is assumed to be constant in this study. Additionally, $Q_{Ai}$ is the flow from the connecting hose to the ith chamber. The cylinder friction is expressed by Equation (19):

$$F_{fric}=\left\{\begin{array}{cc}tanh(a{v}_c)\mid {\tilde{F}}_c\mid (1-{\eta }_c),\hfill & \mathrm{if}{F}_c{v}_c>0\hfill \\ tanh(a{v}_c)\mid {\tilde{F}}_c\mid (\frac{1}{{\eta }_c}-1)\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(19)

where a is the coefficient used to smooth the friction curve versus velocity. ${\eta }_c$ is a constant efficiency of the cylinder.

3.2. The Hoses

The hoses connected between the cylinder and the manifold valves are modeled by Equations (20) and (21) which refer to [19]:

$$\begin{array}{cc}\hfill {\dot{Q}}_{out}& =\frac{(p_1-p_2)A_{hose}-p_f(Q_{out})A_{hose}}{\rho l_{hose}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\dot{p}}_1& =\frac{(Q_{in}-Q_{out})\beta }{A_{hose}{l}_{hose}}\hfill \end{array}$$

(21)

where $Q_{in}$ and $Q_{out}$ are the fluid flows in and out of the hose, and $p_1$ and $p_2$ are the pressures of the inlet and outlet of the hose, respectively. $A_{hose}$ is the area of the hose, $l_{hose}$ is the length of the hose, $\rho$ is the fluid density, and $p_f(Q_{out})$ is the pressure drop across the hose. The pressure drop across a straight pipe/hose can be modeled by the equation:

$$\begin{array}{cc}{p}_{\lambda }\hfill & =\frac{0.3164l_{hose}\rho }{2R{e}^{0.25}{d}_{hose}}\frac{Q_{out}\mid Q_{out}\mid }{{(0.25d_{hose}^2\pi )}^2}(0.5+0.5tanh(\frac{2300-Re}{100}))\hfill \\ & +\frac{128\nu \rho l_{hose}{Q}_{out}}{\pi d_{hose}^4}(0.5+0.5tanh(\frac{-2300+Re}{100}))\hfill \end{array}$$

(22)

where $\nu$ is the kinematic viscosity of the fluid. $Re$ represents the Reynold number which can be computed by Equation (23):

$$Re=\frac{v_{out}{d}_{hose}}{\nu }$$

(23)

Equation (22) combines the pressure loss of the laminar flow and the turbulent flow by the hyperbolic-tangent expression. Consequently, a continuous transition of the pressure loss between the laminar and turbulent flow can be created. When the Reynold number is less than 2200, $(0.5+0.5tanh(\frac{2300-Re}{100}))$ is close to zero, which means the pressure drop is contributed by the laminar flow. On the other hand, when the Reynold number is greater than 2400, $(0.5+0.5tanh(\frac{-2300+Re}{100}))$ is close to zero, which means the pressure drop is contributed by the turbulent flow. Another source of pressure drop is the fitting loss, which can be computed by Equation (24):

$$p_{\zeta }=\zeta \frac{\rho }{2}{Q}_{out}\mid Q_{out}\mid \frac{1}{{(0.25d_{hose}^2\pi )}^2}$$

(24)

where $\zeta$ is the friction coefficient for a given fitting type. Finally, the total resistance in the hose with n line pieces and m fittings can be computed by Equation (25):

$$\begin{array}{cc}\hfill p_f(Q_{out})& =p_{\lambda ,1}(Q_{out})+\dots +p_{\lambda ,n}(Q_{out})+p_{\zeta ,1}(Q_{out})+\dots +p_{\zeta ,m}(Q_{out})\hfill \end{array}$$

(25)

In this paper, the pressure loss of the hoses is modeled by the Equation (26):

$$p_f(Q_{out})=p_{\lambda }(Q_{out})+p_{\zeta ,M}(Q_{out})+p_{\zeta ,C}(Q_{out})$$

(26)

where $p_{\zeta ,M}$ represents the fitting resistance which considers the internal pressure drops in the manifold and $p_{\zeta ,C}$ represents the cylinder inlet loss.

3.3. The Directional Valves

The two-way two-position directional valves are used in this model. The flow across the valve can be described by the following orifice Equation (27):

$$Q_v=sign(\Delta p)C_d{A}_v(\alpha )\sqrt{\frac{2}{\rho }\mid \Delta p\mid }$$

(27)

where $\Delta p$ is the pressure difference cross the valve, $C_d$ is the discharge coefficient, and $A_v(\alpha )$ is the opening area which can be computed by Equations (28)–(30):

$$A_v(\alpha )=\alpha A_0$$

(28)

$$\begin{array}{cc}\hfill \dot{\alpha }& =\left\{\begin{array}{cc}\frac{1}{t_v},\hfill & \mathrm{if}{u}_v=1\hfill \\ -\frac{1}{t_v},\hfill & \mathrm{if}{u}_v=0\hfill \end{array}\right.\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill 0& ⩽\alpha ⩽1\hfill \end{array}$$

(30)

where $A_0$ is the maximum opening area of the valve. In this paper, a total of eight valves are used to control the actuator force. The $t_v$ represents the opening and closing time of the valve. The shifting algorithm of the valve applied in this paper can be further improved [29] to reduce the pressure oscillations and improve the energy efficiency. Moreover, the different opening time has a significant impact on the cylinder pressure [30]. The opening time $t_v$ is selected to be 30 ms in this paper to avoid the cavitation or pressure spikes and to have a relatively fast response to the reference control command. Since the focus of this paper is to examine the controllers’ performance practically, the influence of different opening times is not investigated in this paper.

3.4. The Pressure Accumulators

The accumulators in the DDC system are used as pressure sources and also for energy storage. The dynamics of the pressure accumulator can be modeled with the Equations (31)–(33) [19]:

$$\begin{array}{cc}\hfill {\dot{p}}_{acc}& =\frac{Q_{acc}+\frac{1}{1+\frac{R}{C_v}}\frac{V_g}{T}\frac{1}{{\tau }_a}(T_w-T)}{\frac{V_{a0}-V_g+V_{ext}}{\beta }+\frac{1}{1+\frac{R}{C_v}}\frac{V_g}{p_{acc}}}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\dot{V}}_g& =-Q_{acc}+{\dot{p}}_{acc}\frac{V_{a0}-V_g+V_{ext}}{\beta }\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \dot{T}& =\frac{1}{{\tau }_a}(T_w-T)-\frac{RT}{C_v{V}_g}{\dot{V}}_g\hfill \end{array}$$

(33)

where $p_{acc}$ is the pressure of the accumulator, $Q_{acc}$ is the inlet flow to the accumulator, R is the ideal gas constant, $C_v$ is the gas specific heat at constant volume, $T_w$ is the wall temperature, ${\tau }_a$ is the thermal time constant, $\beta$ is the bulk modulus of the fluid in the pipeline volume $V_{ext}$, $V_{a0}$ is the size of the accumulator, $V_g$ is the gas volume, and T is the gas temperature. Hence the state of the accumulator contains the pressure, the gas volume and the gas temperature. Initially, the state can be specified based on the standard gas law by Equation (34):

$$V_g=\frac{T}{T_0}\frac{p_{a0}}{p_a}{V}_{a0}$$

(34)

where $p_{a0}$ is the pre-charged pressure of the gas at the temperature $T_0$.

3.5. The Hydraulic Motor

For the system presented in this paper, there are 4 chambers and 2 different pressures: the high pressure and the low pressure. The hydraulic motor is connected between the high pressure accumulator and the low pressure accumulator. The flow of the hydraulic motor can be modeled by Equation (35):

$$Q_M=D_w{\omega }_M-\Delta p{C}_{Q1}$$

(35)

where $D_w$ is the displacement of the hydraulic motor, which is constant for a fixed displacement motor, $\Delta p$ is the pressure across the motor, $C_{Q1}$ is the coefficient of the flow loss of the motor, and ${\omega }_M$ is the rotational speed of the motor which is defined by Equation (36):

$${\omega }_M=\frac{p_{avg,exp}\psi }{p_H{k}_{gen}{D}_M}$$

(36)

where $p_{avg,exp}$ is the expected average power output, $p_H$ is the pressure of the high pressure accumulator, $k_{gen}$ is the number of generators, $D_M$ is the total motor displacement, and $\psi$ is a coefficient for the motor speed control to prevent the high pressure from depletion or saturation which is formulated by Equations (37) and (38):

$$\begin{array}{cc}\hfill k& =\frac{4}{(p_{H,max}-p_{H,min})}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \psi & =\left\{\begin{array}{cc}k(p_H-p_{H,min}),\hfill & \mathrm{if}{p}_H>p_{H,min}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(38)

To achieve the desired motor speed introduced in Equation (36), the generator torque control needs to be included. In this paper, the generator and inverter are not modeled and the desired motor speed is assumed achievable. The power in the hydraulic motor can be computed by Equation (39):

$$P_M=\Delta p{Q}_M$$

(39)

This completes the modeling of DDC hydraulic system; the control algorithm is introduced in the next section.

4. The Control Algorithm

Two parts will be presented in this section: the control method for the buoy and the force-shifting algorithm for controlling the valves. The control method for controlling the buoy computes a reference value for the control force at each time step. This reference control force is then used as an input to the PTO , and the actual control force that results from the PTO is computed using the force-shifting algorithm. Each of the two parts is detailed below.

4.1. The Buoy Control Method

Several control methods will be tested in this paper using a simulator that simulates the PTO unit. Some of these controller were originally developed for heave control. It is relatively straightforward, however, to extend a control method from the heave motion to the pitch motion. For example, the singular arc (SA) control method [31] can be used to compute the control torque by the following Equation (40):

$${\tau }_{PTO}(s)=\frac{N(s)}{D(s)}$$

(40)

where:

$$\begin{array}{cc}\hfill N(s)& =(J_{total}{s}^2+(C_r{(sI+A_r)}^{-1}{B}_r-D_r)s+K_{res}){\tau }_{ex}(s)\hfill \\ \hfill D(s)& =s(C_r{(sI+A_r)}^{-1}{B}_r-C_r{(sI-A_r)}^{-1}{B}_r-2D_r)\hfill \end{array}$$

(41)

where the excitation torque can be expressed as Fourier series expansion by Equation (42):

$${\tau }_{ex}=\sum _{i=1}^n\Re ({\tau }_{c,ex}({\omega }_i)\eta ({\omega }_i)e^{i(-{\omega }_{i}t+{\varphi }_i)})$$

(42)

An inverse Laplace transformation is then applied to the SA control to obtain the control in the time domain. The required information to compute the control is the time t, the excitation torque coefficient ${\tau }_{c,ex}$, the wave frequency vector $\overrightarrow{\omega }$ and the time domain phase shift vector $\overrightarrow{\varphi }$.

A reference control method is the feedback proportional-derivative (PD) control. The PD control takes the form of Equation (43):

$${\tau }_{PTO}=K\theta +B\dot{\theta }$$

(43)

where K is the proportional gain and B is the derivative gain.

In addition to the above two control methods, simulated in this paper are model predictive control (MPC) [32], shape-based (SB) control [33], proportional-derivative complex conjugate control (PDC3) [34], and pseudo-spectral (PS) control [35]. Each one of these methods is well documented in the literature, so the details of each control methods are avoided in this paper.

In the original developments, the SA control and the PDC3 control compute a control force that is equivalent to the complex conjugate control (C3), and hence the maximum possible harvested energy in the linear domain. However, the C3 does not account for constraints on the buoy displacement. In fact, since the C3 criterion is to resonate the buoy with the excitation force, the motion of the buoy always violates displacement constraints when controlled using the SA and PDC3 controls. On the other hand, the MPC, SB, and PS control methods compute a control force in an optimal sense, taking displacement constraints into account. Figure 5 shows a simulation for 5 min for the above six control methods when a constraint on the buoy displacement is assumed. The simulation parameters are detailed in Section 5. This simulation does not account for the PTO dynamics and is here presented to highlight the impact of including the PTO into the simulations in Section 5. As can be seen from Figure 5, among the six control methods, the MPC and PD controls performed best, then the SB method, then the PS, and then the PDC3 and SA methods. The two methods (SA and PDC3) that perform best without displacement constraints actually have the poorest performance when accounting for the constraints.

4.2. The Force-Shifting Algorithm

The force-shifting algorithm (FSA) is introduced in this section. The FSA used in this paper is described by Equation (44):

$$\{F_c(t)=\overrightarrow{F}\left[k\right]\mid k=argmin\mid F_{ref}(t)-\overrightarrow{F}\left[k\right]\mid \}$$

(44)

where $F_{ref}$ is the reference control force (computed for instance using one of the six control methods described above), and $\overrightarrow{F}$ is the vector of the possible discrete values for the force. With different permutations of valve openings, it is possible to produce different levels of constant forces, as shown in Figure 6, where it is assumed that $p_H=200$ bar and $p_L=20$ bar. The FSA selects the discrete force level that is closest to the reference control force. It is noted here that the discrete force changes over time due to the fluctuation of the pressures in the accumulators.

5. Simulation Tool

A tool for simulating the dynamics of the WEC including the motion dynamics, the hydrodynamic/hydrostatic force calculations, and the PTO hardware model was developed in MathWorks Simulink${}^{\circledR }$. The detailed Simulink model of the wave energy conversion system is shown in Figure 7. The Plant block simulates the dynamics of the buoy. The PTO block simulates all the equations of the valves, hoses, and accumulators. As can be seen in Figure 7, the excitation force is an input that is computed outside the Plant block. The control force command is computed in the block ‘Control Command’. Despite the name, six different controllers were tested in the ‘Control Command’ block. A detailed Simulink model of the hydraulic system is shown in Figure 8, this model is inside the PTO block in Figure 7.

The parameters of the dynamic model of the wavestar used in the simulations in this paper are listed in Table 1.

5.1. Wave Model

Irregular waves are simulated in this study using the stochastic Pierson–Moskowitz (PM) spectrum. The spectral density is defined by Equation (45):

$$S(\omega )=\frac{5H_{m0}^2}{{\omega }_p{\omega }_n^5}{e}^{-\frac{5}{4{\omega }_n^4}}$$

(45)

where ${\omega }_n=\omega /{\omega }_p$, ${\omega }_p$ is the peak frequency, and $H_{m0}$ is the significant height of the wave. The wave used in the simulation has a significant height of $1.75$ m and a peak period $5.57$ s. The reason to select the PM spectrum is because it is a more conservative choice by considering the power absorption estimation. The JONSWAP spectrum, which is frequently applied, has a narrow frequency band. However, the PM spectrum has a wider frequency band which makes the energy extraction more difficult. Consequently, the PM spectrum is applied in this paper to prevent overestimation of the performance of the controller. With regard to the significant height and peak period applied in the simulation, they are selected from the validated wave climate in [9,36] for the Wavestar C5. The proper range so that the wave for C5 can have major energy absorption is with a significant height from $0.5$ m to 3 m and peak period from 2 s to 7 s. The wave climate applied in this paper is in the middle of the range to avoid losing generality.

5.2. The System Losses

The system losses are computed in this study. The system losses include the pressure loss of the hoses, the flow loss of the generator, and the friction of the cylinder. The pressure loss is shown in Figure 9, in which the vertical line represents the transition between the laminar flow and the turbulent flow when the Reynold number is $Re=2300$, for each of the two possible directions of the fluid flow. The amount of the flow loss and the friction force of the cylinder are shown in Figure 10 and Figure 11. All the system parameters used in the simulations in this paper are listed in

Table 2. The data used in the simulation of the overall WEC system.

Symbol	Value	Unit
Length of the arms
$l_2$	3	m
$l_3$	$2.6$	m
$l_4$	$1.6$	m
Length of the hoses C2M
$l_{A1}$	1	m
$l_{A2}$	1	m
$l_{A3}$	1	m
$l_{A4}$	1	m
Diameter of the hoses C2M
$d_{A1}$	$1.5$	in
$d_{A2}$	$1.5$	in
$d_{A3}$	$1.5$	in
$d_{A4}$	$1.5$	in
Maximum stroke
$x_{c,max}$	3	m
Area of the chambers
$A_1$	$113.4\times {10}^{-4}$	m${}^2$
$A_2$	$32.55\times {10}^{-4}$	m${}^2$
$A_3$	$80.85\times {10}^{-4}$	m${}^2$
$A_4$	$162.75\times {10}^{-4}$	m${}^2$
Max Area of the valves
$A_{01}$	$1.6\times {10}^{-4}$	m${}^2$
$A_{02}$	$1.6\times {10}^{-4}$	m${}^2$
$A_{03}$	$1.6\times {10}^{-4}$	m${}^2$
$A_{04}$	$1.6\times {10}^{-4}$	m${}^2$
Accumulator size
$V_{a0}$	$100\times {10}^{-3}$	m${}^3$
Pressure drop coef
${\zeta }_M$	$1.3$
${\zeta }_C$	1
Specific time constant S
${\tau }_l$	23	s
${\tau }_h$	34	s
Initial pressure of the accumulators
$p_{a,l}$	20	bar
$p_{a,h}$	130	bar
Initial angle
${\alpha }_0$	$1.0821$	rad
Control parameters
K	$-9.16\times {10}^6$	Nm/rad
B	$4.4\times {10}^6$	Nms/rad
Valve opening time
$t_v$	$30\times {10}^{-3}$	s
Wall temperature
$T_w$	50	oC
Ideal gas constant
R	276	J/kg/K
Gas specific heat at constant volume
$C_v$	760	J/kg/K
Motor displacement
$D_w$	100	cc/rev
Flow loss coefficient
$C_{Q1}$	$5.4\times {10}^{-12}$	m${}^3$/s/Pa
Fluid bulk modulus
$\beta$	$1.5\times {10}^9$	Pa

.

6. Simulation Results

The above Simulink tool is used to simulate the performance of the above six control methods. The energy extracted by method is shown in Figure 12. In this simulation, there are limitations on the maximum stroke and the maximum control force. In addition, the PTO dynamics are simulated. The maximum control force in the cylinder in the simulations presented in this paper is assumed to be 215 kN. The maximum allowable displacement in the simulations presented in this paper is assumed to be $1.2$ m. As can be seen in Figure 12, the MPC and PD control methods harvest the highest energy level compared to the other methods. The SB method comes next. The SA, PDC3, and PS control methods come next, and the three of them perform about the same. For further analysis of the performance of the controllers, the capture width ratio (CWR) is evaluated:

$$\begin{array}{c}\hfill CWR=\frac{P_{ave}}{D{P}_w}\end{array}$$

(46)

where $P_{ave}$ is the average power extraction of the buoy, $P_w$ is the wave energy transport, and D is the characteristic dimension of the buoy. The CWRs of the MPC, PD, SB, SA, PDC3, and PS controllers are $51.21\%$, $50.92\%$, $44.68\%$, $37.60\%$, $37.14\%$, and $36.86\%$ respectively. According to the [36], the CWR of the performance of the floater with the applied wave ranges from $40\%$ to $50\%$. Hence, the performance of the proposed controllers is in the reasonable range in terms of energy extraction. Comparing Figure 12 to Figure 5 we can see that by including the PTO model, the performance of the SA method improves slightly, while the performance of the PS degrades slightly, and as a result the three methods PS, SA, and PDC3 perform about the same. The performance of the MPC, PD, and SB control methods actually slightly improve when the PTO model is included.

Another important result to examine is the output mechanical power at the actuator and the output power from the generator. These two quantities are compared in Figure 13. From the figure we can tell that the power absorbed in the generator side is much smoother than the power extracted by the actuator. The hydraulic accumulators act as a power capacitor for energy storage, resulting in this relatively smooth power profile at the generator output. As can be seen in the figure also, the actuator power includes reactive power; these are the times at which the actuator power is negative. At these times, the PTO actually pumps power into the ocean through the actuator. The generator output power does not have any reactive power, confirming that all the reactive power come from the accumulators.

The efficiency of the system is defined by Equation (47):

$${\eta }_{out}=\frac{P_{gen}}{P_{actuator}}$$

(47)

The efficiency depends on the control method. For example, in this test case, the efficiency of the SB controller is $80.15\%$, for the MPC it is $72.58\%$, for the PS it is $67.34\%$, for the SA it is $64.36\%$, and for the PD controller it is $71.76\%$, over 300 s.

In the context of comparing the performance of different control methods, it is important to highlight one significant difference between them that emanates from the theory behind each control method. Each of the MPC, SB, and PS control methods requires wave prediction. That is, wave information (or excitation force) is needed over a future horizon at each time step in the simulation. In the simulations in this paper, this future horizon is assumed to be $0.6$ s for the SB and MPC control methods, and is assumed to be 60 s for the PS control. Wave prediction is assumed perfect in these simulations. Non-perfect wave prediction would affect the results obtained using these methods. The PD, SA, and PDC3 control methods do not need future wave prediction.

This simulation tool also provides detailed operation information that is useful for characterizing different components in the system. For example, the generator speed is computed in the simulation, and is shown in Figure 14. As shown in the figure, the speed is oscillating around 1200 RPM.

To present detailed plots for the response of the buoy, only one control method is selected as a sample to avoid excessive figures in the paper. The SB method is selected here to present the detailed WEC response in this section. The angular displacement of the buoy is shown in Figure 15; the maximum angular displacement is about 10 degrees and it is below 5 degrees most of the time. The angular velocity of the buoy is shown in Figure 16. The cylinder force and the PTO torque are shown in Figure 17 and Figure 18, respectively. Both the reference and actual values are plotted in each of the two figures. As can be seen in Figure 17, the control force is below the force limit of 215 kN. The accumulator pressure is shown in Figure 19. The high pressure is oscillating around 100 bar, while the low pressure is stable around 20 bar. The chamber pressure is shown in Figure 20. Significant fluctuations can be observed when the hydraulic system is extracting energy. This is necessary to be able to track the reference control command effectively. However, those fluctuations may be reduced by increasing the valve opening area or including relief valves.

7. Discussion

In this paper, different recent control methods are tested using a simulation tool that simulates a hydraulic PTO system. In a theoretical test (where PTO is assumed to track the reference control command ideally and in the absence of all constraints,) the SA controller has the best performance in terms of energy extraction. However, the performance of the SA controller with the hydraulic system model included is the poorest among the tested six control methods. To get more insight into this phenomenon, consider

Table 3. Capacity requirement of the controllers without hydraulic system.

Symbol	Value	Unit
The SA controller
$F_{PTO,max}$	3705	kN
$x_{c,max}-x_{c,min}$	$3.2$	m
$P_{ave}$	$35.11$	kW
The PD controller
$F_{PTO,max}$	1119	kN
$x_{c,max}-x_{c,min}$	$1.1$	m
$P_{ave}$	$21.00$	kW
The PDC3 controller
$F_{PTO,max}$	1404	kN
$x_{c,max}-x_{c,min}$	$1.6$	m
$P_{ave}$	$21.08$	kW

that presents data for three controllers (SA, PD, and PDC3) in the theoretical test case. As can be seen in Table 3, the energy extracted by the PD controller in this theoretical test is about $60\%$ of that of the SA controller. However, the buoy maximum displacement associated with the SA control is significantly greater than that of the PD control (almost three times higher) which makes it more difficult to achieve. Similarly, the maximum control force required by the SA control is significantly greater than that of the PD control, which means a PTO might not be able to track the command force at all times when using a SA control, while it is more likely to track a command force generated using a PD control. The data of the PDC3 control in Table 3 also highlights that the PDC3 control in this test case generates about the same level of average power, but in a higher displacement range and with higher force capability. This indicates that including a model for the PTO would result in favorable performance for the PD control compared to the PDC3. To highlight the impact of the PTO model on the performance of the different control strategies, consider

Table 4. Capacity requirement of the controllers with hydraulic system.

Symbol	Value	Unit
The SA controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$0.96$	m
$P_{ave}$	$13.49$	kW
The PD controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$1.1$	m
$P_{ave}$	$18.26$	kW
The PDC3 controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$1.2$	m
$P_{ave}$	$13.32$	kW
The SB controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$0.8$	m
$P_{ave}$	$16.02$	kW
The MPC controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$1.1$	m
$P_{ave}$	$18.37$	kW
The PS controller
$F_{PTO,max}$	215	kN
$x_{c,max}-x_{c,min}$	$0.90$	m
$P_{ave}$	$13.22$	kW

. The data are presented for all the six control methods. As can be seen from Table 4, all the control methods reached the maximum possible control capacity allowable by the PTO. Since this maximum control force is well below that needed by the SA in Table 3, the amount of harvested energy in this practical case is significantly less than that computed in the theoretical case ($13.49$ kW compared to $35.11$ kW in average power). The drop in energy harvested using the PD control however is less since the maximum force needed theoretically was as high as that of the SA. The displacement of the PDC3 reached the maximum displacement allowable by the WEC ($1.2$ m.) This is expected since the PDC3 tends to increase the displacement and hence it would reach a limit imposed by the WEC system.

8. Conclusions

The main conclusion of this paper is that a controller that is optimal in theoretical analysis might not be optimal when tested in a practical test environment. In particular, this paper sheds light on considerations that need to be accounted for in designing a control method for WEC systems. The first of these considerations is the limitation of the maximum possible PTO control force. This limitation impacts methods such as the singular arc control which is a control method developed in an optimal sense using classical optimal control theory. The second consideration is the limitation due to the maximum possible displacement of the WEC system. This limitation impacts the optimality of some control methods such as the multi-resonant proportional derivative control that is derived in an optimal sense to satisfy the complex conjugate criterion. Another consideration is the capability of the PTO to track the control command. The hydraulic PTO presented in this paper produces discrete levels of control forces and hence the dynamics of this PTO need to be accounted for in designing a control system for practical energy harvesting.