Hybrid Model Predictive Control of a Two-Body Wave Energy Converter with Mechanically Driven Power Take-Off

Abstract

In this paper, a variable damper is proposed to regulate the efficiency of a two-body wave energy converter (WEC) with mechanically driven power take-off (PTO). The variable damper introduces logic constraints into the WEC system, which can be translated into a mixed logical dynamical form with the dynamics of real-valued variables, the dynamics of logic variables, and their interactions. A hybrid model predictive control (MPC) method is used to determine the control inputs, which has the capacity to handle various constraints. The performance is assessed through simulations to evaluate the effectiveness of the proposed method. The achievable performance improvements of the proposed hybrid MPC are shown by means of comparative analysis with uncontrolled WEC devices. The results show that the proposed hybrid MPC has a high requirement on the lower bound of the variable damper and the maximum damping is used only at low relative velocities to achieve the optimum phase, like latching control. The hybrid MPC performs exceptionally well under wave conditions with a small significant wave height and long wave period, improving the power generation of the uncontrolled system up to 22.5%. And, the prediction error has a significant effect on hybrid MPC performance, especially for long prediction horizon.

1. Introduction

Among renewable energy resources, ocean wave energy has the advantages of high energy density and continuous power supply [1,2]. The tremendous potential of wave energy has attracted much research interest to wave energy conversion technologies. While the breadth of knowledge and technologies has become quite large, the commercial exploitation of wave energy is still non-operational due to the high levelized cost of energy (LCOE).

The implementation of an advanced control strategy is an efficient way to improve conversion efficiency and lower the LCOE of WECs. An ideal control technique is the reactive control method, also known as complex conjugate control, which calls for the system to function as a motor in order to achieve the optimum phase and amplitude of WEC oscillation velocity [3]. It is desirable that the PTO machinery has high energy conversion efficiency to return some energy during part of the oscillation cycle [4]. However, returning the necessary amount of energy is difficult for PTO machinery, which means that reactive control cannot be implemented in practice. Budal and Falnes [5] proposed a latching control method, which can only be applied to devices with a resonant period shorter than the dominant wave period [6]. The latching control is implemented by locking the PTO for certain time intervals to lower the resonant frequency of the device to match the incident wave frequency. Declutching control is a counterpart to latching control, which is suitable for devices with a resonant period longer than the dominant wave period [7]. The declutching control method is implemented by bypassing the PTO mechanism for certain time intervals to increase the resonant frequency of the device to match the incident wave frequency. Both latching control and declutching control are not optimal control strategies for WEC devices since they simply seek to reach the ideal phase conditions while ignoring the optimum amplitude. A combination of latching control and declutching control that can switch between power generation, latching mode, and declutching mode was proposed by Feng and Kerrigan [8]. The latching–declutching control can extend the range of wave periods to make the WEC generate more electricity. These control methods are easily implemented in idealized cases for monochromatic waves and inactive physical constraints. However, they are too complicated to be implemented in realistic scenarios because the device performance is difficult to be evaluated in terms of the requirements of the method [9] and the latching or declutching time is difficult to be determined for various sea states.

The development of advanced control strategies for constrained WECs is necessary to improve the energy harvest as much as possible while maintaining safe operation. The energy maximization for WEC devices can be regarded as a constrained optimal control problem, which is noncausal [10] and needs to be combined with a wave excitation force prediction algorithm. MPC is suitable for the WEC control problem since various constraints can be easily handled, and the prediction of wave excitation can be explicitly incorporated. Lin and Huang [11] proposed a reactive rollout MPC formulation that introduces a terminal value function and a terminal constraint set to improve the energy extraction of WEC by taking into account additional future wave information. The proposed method has significant superiority over the traditional MPC, especially in the case of a short optimization horizon. To improve the accuracy of the WEC model and reduce the computational burden, Li et al. [9] proposed an improved MPC with sliding mode control, where the sliding mode control is employed to compensate for the model mismatch. The improved MPC is widely feasible for different processers and sea conditions by adjusting the proper parameters while maintaining high energy conversion efficiency. Sergiienko et al. investigated [12] the effect of MPC parameters, including prediction horizon and control force penalties, on the design of a PTO system. Approximate control force penalty terms can significantly reduce the requirements for peak control force and slew rate. Faedo et al. [13] provided a critical comparison of the MPC and MPC-like algorithms, while the MPC algorithms within the broader context of other “optimal” control schemes were presented.

In the previously mentioned investigations, it is assumed that the PTO force is controllable, and the PTO force is usually regarded as the optimization variable in the optimal control problem of WEC. However, there are differences in how accurately forces can be achieved for the wide range of PTO systems. The hydraulic PTO system is frequently employed in the design of WEC due to its advantages of having high technology maturity and a huge power load capacity [14]. However, hydraulic leakage, which is the primary concern for the marine environment, and the costly regular maintenance as a result of the complexity of hydraulic transmission are the drawbacks of the hydraulic PTO system [15]. Another type of PTO system is the pneumatic air turbine, mainly used in oscillating water column devices. The fact that key components of the air turbine PTO system are situated above water significantly enhances the reliability and maintainability of the WEC device. However, the main drawbacks come from the challenges associated with the compressibility of the air and the directional airflow, which make the system less efficient in energy conversion [16]. The hydraulic motor and pneumatic air turbine transmission systems require a fluid power transfer medium to accomplish the energy conversion, which is inevitably accompanied by energy loss. The direct mechanical drive system converts the mechanical energy directly into electricity with up to 97% efficiency [17]. Li et al. [18] designed a self-reactive WEC based on a ballscrew PTO. Martin et al. [19] analyzed and validated a two-body WEC with a ballscrew using numerical methods and wave tank testing. A key component of a mechanically driven PTO is the mechanical motion rectifier, which can convert bidirectional wave motion into a unidirectional input at the generator shaft. However, one main drawback of the mechanical drive PTO is that it is difficult to be controlled. Yang et al. [20] proposed a new PTO based on an active mechanical motion rectifier, which can achieve the controllability of the WEC device while guaranteeing unidirectional transmission. Li et al. [21] proposed a variable inertia flywheel to directly manipulate the equivalent mass of the WEC device with the aim of altering the system resonance frequency and power absorption bandwidth. The numerical studies showed that the variable inertia wheel has the capacity to adapt to various wave conditions and increase power output. Yang et al. [22] investigated the meshing and overrunning phases due to the introduction of a variable inertia flywheel with a passive configuration and verified through experiments that the variable inertia wheel is beneficial in improving energy conversion efficiency.

In this study, the realistic feasibility of the control of the WEC is considered and a hybrid MPC method is proposed. The control mechanism is a variable damper, which has proven applications in various industries and is placed in parallel with the PTO in the WEC system. In order to verify the superiority of the proposed method, the hybrid MPC method and no-control WEC are simulated and analyzed under different wave conditions. Additionally, the control performances under various parameter conditions are compared to explore the influence and function of parameters in the proposed strategy on the two-body WEC.

The rest of the paper is organized as follows. In Section 2, the forces on the two-body WEC are analyzed, and the state-space model is established. Section 3 introduces the proposed hybrid MPC method. Section 4 investigates the effects of constraints and parameters in the proposed method on the two-body WEC. Section 5 concludes the research.

2. Wave Energy Converter Model

2.1. Structure

The two-body point absorber WEC considered in this paper is shown in Figure 1. It consists of a float, a reaction section, a PTO system, and a control system. The reaction section consists of a spar and a heave plate. The PTO system consists of a ballscrew, a ballscrew nut, a gearbox, and a generator. The relative motion between the float and reaction section, driven by incident waves, is converted into rotation motion of the PTO system through the inverse driven mode of the ballscrew and ballscrew nut. The rotation motion is amplified by the gearbox and drives a generator to generate electricity into the battery. The control system changes the internal damping of the WEC through a variable damper, which is placed in parallel with the PTO system.

Only the heave motion of WEC is considered in this study. The equations of motion for the WEC can be expressed as follows [23]:

$$\begin{array}{l}{m}_1{\ddot{x}}_1=F_{exc1}+F_{rad11}+F_{rad12}+F_{b1}+F_{pto}+F_u\\ m_2{\ddot{x}}_2=F_{exc2}+F_{rad21}+F_{rad22}+F_{b2}-F_{pto}-F_u\end{array}$$

(1)

where subscripts “1” and “2” represent the float and reaction section, respectively; m is the dry mass; x is the heave displacement with respect to the equilibrium position; F_exc represents the wave excitation force; F_radij is the radiation force imposed on the body i by the oscillation of the body j; F_b is the net buoyancy restoring force; F_pto is the PTO force; and F_u is the control force. The wave excitation force F_exc, radiation force F_rad, and net buoyancy restoring force F_b are hydrodynamic forces, which can be calculated using hydrodynamic coefficients provided by frequency-domain boundary element method (BEM) solvers, such as WAMIT, AQWA, and NEMOH.

2.2. Hydrodynamic Force

The regular wave excitation force can be represented as follows:

$$F_{exc}(t)=\Re \left[\frac{H}{2}{F}_{exc}(\omega )e^{i\omega t}\right]$$

(2)

where $\Re$ denotes the real part of the formula; H is the wave height; and $F_{exc}(\omega )$ is the frequency-dependent complex amplitude of the wave excitation force. The dependence of $F_{exc}(\omega )$ on the wave direction is neglected, considering the symmetry of the WEC device.

The irregular wave excitation force can be calculated as the superposition of regular wave excitation across all wave frequencies, as follows:

$$F_{exc}(t)=\Re \left[{\displaystyle \sum _{j=1}^N{A}_j}{F}_{exc}\left({\omega }_j\right)e^{i\left({\omega }_{j}t+{\varphi }_j\right)}\right]$$

(3)

where A_j is the wave amplitude at angular frequency ${\omega }_j$, which can be obtained from the wave spectrum; ${\varphi }_j$ is the randomized phase angle at angular frequency ${\omega }_j$, and N is the number of frequency bands selected to discretize the wave spectrum. The Pierson–Moskowitz (PM) spectrum was chosen to describe incident waves in this paper.

$$S_{PM}(f)=\frac{H_s^2}{4}{\left(1.057f_p\right)}^4 f^{-5}\exp \left[-\frac{5}{4}{\left(\frac{f_p}{f}\right)}^4\right]$$

(4)

$$A_j=\frac{H_j}{2}=\sqrt{2S\left({\omega }_j\right)d{\omega }_j}$$

(5)

where f is the frequency of the PM spectrum; H_s is the significant wave height; f_p is the peak frequency; and dω_j is the frequency interval.

The radiation force can be represented as follows [23]:

$$F_{radij}(t)=-{\mu }_{ij}{\ddot{x}}_j-{\displaystyle {\int }_0^t{K}_{rij}}(t-\tau ){\dot{x}}_j(\tau )d\tau$$

(6)

where μ_ij is the added mass at infinite frequency and K_rij is the radiation impulse response function. To improve the computing efficiency, a state-space approximation of the convolution term in Equation (6) is needed. It is constructed as follows:

$$\begin{array}{r}\dot{q}=A_{r}q+B_r\dot{x}\\ {\displaystyle {\int }_0^t{K}_r}(t-\tau )\dot{x}d\tau =C_{r}q+D_r\dot{x}\end{array}$$

(7)

where A_r is the system matrix, the dimension of which depends on the accuracy of approximation representation; B_r is the input matrix; C_r is the output matrix; D_r is the feedthrough matrix, which is assumed to be zero to maintain the causality of the system and will be omitted in the following equations; and q and $\dot{q}$ are the system state and its derivative, respectively. The system input is the heave speed of the corresponding body. Figure 2 shows the impulse response functions and their state-space realizations. The approximate accuracies of all state-space realizations are above 95%. The dimensions of all system matrices are $2\times 2$.

The net buoyancy restoring force can be represented as follows:

$$F_b=-\rho gSx=-K_{H}x$$

(8)

where ρ is the seawater density; g is the gravity acceleration; S is the water plane area of the body; and K_H is the hydrostatic stiffness coefficient.

2.3. Power Take-Off Force

The transmission of the PTO system consists of three main parts, the translation of linear motion to rotational motion, rotation speed amplification, and electricity generation. The dynamics of the ballscrew can be considered as follows:

$$\omega =2\pi \frac{{\dot{x}}_1-{\dot{x}}_2}{L}$$

(9)

where ω is the ballscrew angular speed, that is, the input angular speed of the gearbox; ${\dot{x}}_1-{\dot{x}}_2$ is the relative heave speed between the float and reaction section; and L is the ballscrew lead. The relation between PTO force and the input moment of the gearbox can be described as follows:

$$M_i=-\frac{F_{pto}}{2\pi }L$$

(10)

where M_i is the input moment of the gearbox, and the transmission ratio of the gearbox is denoted as n_g. The angular speed and input moment of the generator are calculated as follows:

$$\begin{array}{c}{\omega }_g=n_g\omega \\ M_g=\frac{1}{n_g}{M}_i\end{array}$$

(11)

where ω_g is the angular speed of the generator, and M_g is the input shaft moment of the generator, which is the sum of the electric torque load and the torque induced by the inertia of the generator.

$$M_g=M_r+I_g{\dot{\omega }}_g$$

(12)

where M_r is the electric moment induced by the current in the energy harvesting circuit, and I_g is the moment of inertia of the generator.

$$M_r=k_{t}i=\frac{k_{t}V}{R_i+R_o}=\frac{k_t{k}_e}{R_i+R_o}{\omega }_g$$

(13)

where k_t is the torque constant of the generator; i is the output current of the generator; V is the output voltage of the generator; R_i is the internal resistance (motor armature resistance) in the circuit; and R_o is the external resistance (passive load) in the circuit [24,25,26].

Substituting Equations (11)–(13) into Equation (10), the PTO force can be calculated as follows.

$$n_g^2{I}_g\frac{4{\pi }^2}{L^2}({\ddot{x}}_1-{\ddot{x}}_2)+\frac{4{\pi }^2{n}_g^2{k}_t{k}_e}{L^2(R_i+R_o)}({\dot{x}}_1-{\dot{x}}_2)=-F_{pto}$$

(14)

The PTO mechanism can be represented as a linear mass-spring-damper system, where the PTO force is given by the following equation.

$$F_{pto}=-m_{pto}({\ddot{x}}_1-{\ddot{x}}_2)-c_{pto}({\dot{x}}_1-{\dot{x}}_2)-k_{pto}(x_1-x_2)$$

(15)

where m_pto is the mass of the PTO system; c_pto is the damping of the PTO system; and k_pto is the stiffness of the PTO system. By putting together Equations (14) and (15), it can be deduced as follows.

$$\{\begin{array}{l}{m}_{pto}=n_g^2{I}_g\frac{4{\pi }^2}{L^2}\\ c_{pto}=\frac{4{\pi }^2{n}_g^2{k}_t{k}_e}{L^2(R_i+R_o)}\\ k_{pto}=0\end{array}$$

(16)

2.4. Control Force

In the proposed controller, the control force changes with the variable damping, and the damper force can be represented as follows:

$$F_u=\beta ({\dot{x}}_1-{\dot{x}}_2)$$

(17)

where β is the controllable damping coefficient. Variable damping can be obtained via an electrohydraulic damper, whereby the concept of which is depicted in Figure 3.

Compared to the classical passive hydraulic damper, the electrohydraulic damper comprises electronic valves instead of passive valves. The required damping coefficient can be achieved by adjusting the opening aperture of electronic valves appropriately, controlled by the current i of the servo mechanisms. The energy consumed in obtaining an appropriate damping coefficient is relatively small or negligible. A typical force–current map for an electrohydraulic damper is shown in Figure 4. It should be noted that the positive direction of damper force F_d is opposite to that of the control force $F_u$, i.e., F_d = −F_u.

It can be seen that the damper force non-linearly depends on the relative velocity between rod, cylinder, and current, and is constrained within a shadow-bounded region consisting of five straight lines.

$$\{\begin{array}{l}u={\beta }_{1}v+{\alpha }_1\\ u={\beta }_{2}v\\ u={\beta }_{3}v\\ u={\beta }_{4}v\\ u={\beta }_{5}v+{\alpha }_5\end{array}$$

(18)

The required control current for driving the variable damper system is based on calculated control force F_u and relative velocity ${\dot{x}}_1-{\dot{x}}_2$ at each time, which can be described as a non-linear inverse function of the damper force, as follows:

$$\{\begin{array}{l}{F}_d=f(i,{\dot{x}}_1-{\dot{x}}_2)\\ i=f^{-1}(F_d,{\dot{x}}_1-{\dot{x}}_2)\end{array}$$

(19)

The detailed function can be obtained experimentally based on a numerical model, such as a polynomial model [27]. In this study, we focus more on the design of the control strategy, and so the primitive function and inverse function are assumed to be known.

2.5. State-Space Model

To control the WEC device via advanced control strategies, a discrete state-space representation for the WEC system is necessary. The state-space formulation of the convolution integral term in radiation force representation for the two-body WEC system was constructed, which can be formulated as follows.

$$\begin{array}{l}\underset{\dot{q}}{\underbrace{\left[\begin{array}{c}{\dot{q}}_{11}\\ {\dot{q}}_{12}\\ {\dot{q}}_{21}\\ {\dot{q}}_{22}\end{array}\right]}}=\underset{A_r}{\underbrace{\left[\begin{array}{cccc}{A}_{r11}& 0& 0& 0\\ 0& A_{r21}& 0& 0\\ 0& 0& A_{r21}& 0\\ 0& 0& 0& A_{r22}\end{array}\right]}}\underset{q}{\underbrace{\left[\begin{array}{c}{q}_{11}\\ q_{12}\\ q_{21}\\ q_{22}\end{array}\right]}}+\underset{B_r}{\underbrace{\left[\begin{array}{cc}{B}_{r11}& 0\\ 0& B_{r12}\\ B_{r21}& 0\\ 0& B_{r22}\end{array}\right]}}\underset{\dot{x}}{\underbrace{\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]}}\\ \left[\begin{array}{c}{F}_{rad1}^c\\ F_{rad2}^c\end{array}\right]=\underset{C_r}{\underbrace{\left[\begin{array}{cccc}{C}_{r11}& C_{r12}& 0& 0\\ 0& 0& C_{r21}& C_{r22}\end{array}\right]}}q\end{array}$$

(20)

where $F_{rad1}^c$ and $F_{rad2}^c$ are the radiation forces (convolution integral term) loading in the float and reaction section, respectively. $q_{\mathit{i}\mathit{j}}$, $i,j\in \left\{1,2\right\}$ is the system state of the state-space realization for the convolution term of radiation forces $F_{radij}$.

$$\begin{array}{l}{F}_{rad1}^c={\displaystyle {\int }_0^t{K}_{r11}}(t-\tau )\dot{x}d\tau +{\displaystyle {\int }_0^t{K}_{r12}}(t-\tau )\dot{x}d\tau \\ F_{rad2}^c={\displaystyle {\int }_0^t{K}_{r21}}(t-\tau )\dot{x}d\tau +{\displaystyle {\int }_0^t{K}_{r22}}(t-\tau )\dot{x}d\tau \end{array}$$

(21)

Substituting Equations (6), (8), (15) and (21) into Equation (1), the dynamic equation can be restated as follows:

$$M\ddot{x}=-\left[\begin{array}{l}{F}_{rad1}^c\\ F_{rad2}^c\end{array}\right]-c_{pto}N\dot{x}-\left(K_H+k_{pto}N\right)x+\underset{F_{exc}}{\underbrace{\left[\begin{array}{c}{F}_{exc1}\\ F_{exc2}\end{array}\right]}}+\left[\begin{array}{c}{F}_u\\ -F_u\end{array}\right]$$

(22)

where,

$$\begin{array}{l}M=\left[\begin{array}{cc}{\mu }_{11}& {\mu }_{12}\\ {\mu }_{21}& {\mu }_{22}\end{array}\right]+\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]+m_{pto}N\\ K_H=\left[\begin{array}{cc}{K}_{H1}& 0\\ 0& K_{H2}\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} N=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}$$

(23)

It can be assumed that the state vector is chosen as $p={\left[\begin{array}{ccccc}{x}_1& x_2& {\dot{x}}_1& {\dot{x}}_2& q^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$; then, the continuous state-space model can be obtained by substituting Equation (20) into Equation (22).

$$\begin{array}{l}\dot{p}=Ap+B_e{F}_{exc}+B_u{F}_u\\ A=\left[\begin{array}{ccc}0& I& 0\\ -M^{-1}(c_{pto}N+K_H)& -M^{-1}(c_{pto}N)& -M^{-1}{C}_r\\ 0& B_r& A_r\end{array}\right]\\ B_e=\left[\begin{array}{c}0\\ M^{-1}\\ 0\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} B_u=\left[\begin{array}{c}1\\ -1\end{array}\right]B_e\end{array}$$

(24)

Finally, the state-space model is discretized by zero-order hold (ZOH) with sampling time T_s (T_s = 0.25 s in this study) to obtain the discrete state-space representation for the WEC system.

$$\begin{array}{c}p(k+1)=A_{d}p(k)+B_{de}{F}_{exc}(k)+B_{du}{F}_u(k)\\ =A_{d}p(k)+B_{de}{F}_{exc}(k)-B_{du}{F}_d(k)\end{array}$$

(25)

where,

$$\begin{array}{l}{A}_d=e^{A{T}_s}\\ B_{de}={\displaystyle {\int }_0^{T_s}{e}^{A{T}_s}{B}_{e}dt}\\ B_{du}={\displaystyle {\int }_0^{T_s}{e}^{A{T}_s}{B}_{u}dt}\end{array}$$

(26)

3. Hybrid MPC Problem Formulation

3.1. Hybrid Dynamical Model

As mentioned in Section 2.4, the damper force must lay in the shadow-bounded region in the force–current map. This non-linear constraint can be translated into a set of thresholds and logic conditions by introducing a binary variable δ_v.

$$\left[{\delta }_v=1\right]\leftrightarrow [{\dot{x}}_1-{\dot{x}}_2\ge 0]$$

(27)

where the connective “$\leftrightarrow$” means “if and only if”. This mixed logic continuous relation can be translated into linear constraints involving the logical variable δ_v.

$$\begin{array}{l}{\dot{x}}_1-{\dot{x}}_2\ge v_l(1-{\delta }_v)\\ {\dot{x}}_1-{\dot{x}}_2\le \epsilon +(v_u-\epsilon ){\delta }_v\end{array}$$

(28)

where v_l and v_u are lower and upper bounds to the function ${\dot{x}}_1-{\dot{x}}_2$, respectively; $\epsilon >0$ is a small tolerance, typically the machine precision.

The constraints for damper force can be written in a more detailed form, as follows [28]:

$$\begin{array}{l}{z}_1=\{\begin{array}{l}{\beta }_2({\dot{x}}_1-{\dot{x}}_2)+{\alpha }_2-F_d,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{if} {\dot{x}}_1-{\dot{x}}_2\ge 0\\ F_d-{\beta }_2({\dot{x}}_1-{\dot{x}}_2)-{\alpha }_2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{otherwise}\end{array}\\ z_2=\{\begin{array}{l}{\beta }_1({\dot{x}}_1-{\dot{x}}_2)+{\alpha }_1-F_d,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{if} {\dot{x}}_1-{\dot{x}}_2\ge 0\\ F_d-{\beta }_5({\dot{x}}_1-{\dot{x}}_2)-{\alpha }_5,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{otherwise}\end{array}\\ z_3=\{\begin{array}{l}{F}_d-{\beta }_4({\dot{x}}_1-{\dot{x}}_2)-{\alpha }_4,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{if} {\dot{x}}_1-{\dot{x}}_2\ge 0\\ {\beta }_3({\dot{x}}_1-{\dot{x}}_2)+{\alpha }_3-F_d,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{otherwise}\end{array}\end{array}$$

(29)

where z₁, z₂, and z₃ $\in ℝ$ are continuous auxiliary variables, which are imposed constraints.

$$z_1,\hspace{0.17em}{z}_2, z_3\ge 0$$

(30)

A generic if–then–else construct of the form is as follows:

$$\begin{array}{l}\mathrm{IF}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \delta \hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{THEN}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} z=a_{1}x+b_{1}u+f_1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{ELSE}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} z=a_{2}x+b_{2}u+f_2\end{array}$$

(31)

where $\delta \in \left\{0,1\right\}$, $z\in ℝ$, $x\in {ℝ}^n$, $u\in {ℝ}^m$, a₁, a₂, b₁, b₂, f₁, and f₂ are constants with approximate dimensions, which can be translated into integer linear inequalities of the form.

$$\begin{array}{l}\left(f_2-F_1\right)\delta +z\le a_{2}x+b_{2}u+c_2\\ \left(f_1-F_2\right)\delta -z\le -a_{2}x-b_{2}u-c_2\\ \left(f_1-F_2\right)(1-\delta )+z\le a_{1}x+b_{1}u+c_1\\ \left(f_2-F_1\right)(1-\delta )-z \le -a_{1}x-b_{1}u-c_1\end{array}$$

(32)

where f_i and F_i are lower and upper bounds to $a_{i}x+b_{i}u+c_i$, $i\in \left\{1,2\right\}$. Each if–then–else equation in Equation (29) can be translated into four linear inequality constraints.

Relative displacement and speed constraints are other issues that need to be considered for WEC devices.

$$\begin{array}{l}{x}_{min}\le x_1-x_2\le x_{max}\\ {\dot{x}}_{min}\le {\dot{x}}_1-{\dot{x}}_2\le {\dot{x}}_{max}\end{array}$$

(33)

where x_min and x_max are the maximum and minimum relative displacement that the system can achieve. ${\dot{x}}_{min}$ and ${\dot{x}}_{max}$ are the maximum and minimum relative velocity. The relative displacement constraint represents a finite stroke limit of the WEC device, which prevents the PTO hardware and the surrounding structure from impulsive loads [29].

The discrete-time dynamic system (25) with constraints (28), (30), (32) and (33) can be modeled as a mixed logical dynamical (MLD) system in HYSDEL (a high-level modeling language for the specification of hybrid systems) [30]. It should be noted that the relative velocity constraints are only used to guarantee the operation of HYSDEL.

$$\begin{array}{l}p(k+1)=A_{d}p(k)-B_{du}{F}_d(k)+B_{d\delta }{\delta }_v(k)+B_{dz}z(k)+B_{de}{F}_{exc}\\ E_2{\delta }_v(k)+E_{3}z(k)\le E_1{F}_d(k)+E_{4}p(k)+E_5\end{array}$$

(34)

In this paper, $p\in {ℝ}^{12}$ is the state vector; $F_d\in ℝ$ is the control input;${\delta }_v\in \left\{0,1\right\}$ is the logic variable; $z\in {ℝ}^3$ is the vector of continuous auxiliary variables; $F_{exc}$ is the wave excitation force vector, which can be seen as the disturbances of the state update equation; A_d, B_du, and B_de are the discrete system matrices from the original state-space Equation (25); B_dδ and B_dz are the zero matrices in our case; and E₂, E₃, E₁, E₄, and E₅ are matrices of suitable dimensions.

3.2. Hybrid MPC Problem

The objective of control in WEC is to maximize the power generation over the prediction horizon, which can be described as a receding horizon optimal problem. The objective function can be written as follows:

$$\underset{u}{\max }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} {\displaystyle \sum _{i=k}^{k+N_p}{c}_{pto}{\left[{\dot{x}}_1(i)-{\dot{x}}_2(i)\right]}^2}$$

(35)

For convenience, the objective function can be written as a cost function (and thus is a minimization problem) of the form.

$$\underset{u}{\min }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} {\displaystyle \sum _{i=k}^{k+N_p}-c_{pto}{p}^{\mathrm{T}}(i)Qp(i)}$$

(36)

where N_p is the optimization horizon length.

$$\begin{array}{l}Q=N^{\mathrm{T}}N\\ N=\left[\begin{array}{ccccc}0& 0& 1& -1& 0^{\mathrm{T}}\end{array}\right]\end{array}$$

(37)

Based on the hybrid logic constraints, the hybrid MPC can be formulated as follows:

$$\begin{array}{l}\underset{u}{\min }{\displaystyle \sum _{i=k}^{k+N_p}{c}_{pto}{p}^{\mathrm{T}}(i)Qp(i)}\\ \mathrm{s}.\mathrm{t}.\hspace{0.17em}\hspace{0.17em}\{\begin{array}{l}p(k+1)=A_{d}p(k)-B_{du}{F}_d(k)+B_{d\delta }{\delta }_v(k)+B_{dz}z(k)+B_{de}{F}_{exc}\\ E_2{\delta }_v(k)+E_{3}z(k)\le E_1{F}_d(k)+E_{4}p(k)+E_5\end{array}\end{array}$$

(38)

4. Numerical Investigation and Result

4.1. Case Study

In this section, a case of hybrid MPC for the two-body WEC is simulated. The key physical specifications for the WEC system are given in

Table 1. Specifications of the WEC model (CoG: center of gravity).

Body		Diameter (m)	Draft (m)	Mass (t)	CoG (m)
Float		20	3	727.01	−0.72
Reaction section	Spar	6	30	878.30	−21.29
	Plate	30

(the parameters are from Reference Model 3 [31]). The seawater density ρ takes the value of 1025 kg/m³ and the gravity acceleration takes the value of 9.81 m/s².

Under the assumption that the future wave excitation force is known, the hybrid MPC was applied in the two-body WEC with T_p = 8 s and H_s = 2.5 m wave input. The spectrum and corresponding wave elevation are shown in Figure 5. The prediction horizon N_p = 32 (8 s). The entire simulation time is 120 s.

The same PTO mechanism was used for all simulations. The parameters of the PTO system can be found in

Table 2. Power take-off parameters.

Symbol	Quantity	Parameter Explanation
ng	10	gearbox ratio
Ig	0.54 kg∙m2	generator moment of inertia
L	0.1 m	ballscrew lead
kt	2.5 Nm/A	generator torque constant
ke	120 Vs/rad	generator speed constant
Ri	2.1 Ω	internal resistance
Ro	100 Ω	external resistance

. The boundary coefficients of damper force are shown in

Table 3. Damper force constraint parameters.

Parameter	β1	β2	β3	β4	β5	α1	α5
Value	1.2 × 106	2.4 × 106	1.0 × 104	1.0 × 104	1.8 × 106	6.0 × 105	−4.2 × 105

, and the relative displacement and velocity constraints are summarized in

Table 4. Relative displacement and velocity constraints.

Symbol	Quantity	Parameter Explanation
xmax	8 m	maximum relative displacement
xmin	−8 m	minimum relative displacement
vmin	5 m/s	maximum relative velocity
vmax	−5 m/s	minimum relative velocity

.

The hydrodynamic coefficients used in this study were provided by AQWA 2022. The wave excitation forces, radiation forces, and buoyancy restoring forces under no control and hybrid MPC are shown in Figure 6 and Figure 7, respectively. It should be noted that the wave excitation forces are the same under uncontrolled and hybrid MPC, due to the limitations of the BEM software.

Figure 8 shows the time history curves of the relative displacement and velocity of the WEC. It can be seen that the motion response of the two-body WEC is dramatically oscillatory. The result of the oscillatory period will be neglected in the following sections. Moreover, we can see that the relative displacement and velocity oscillate around zero and do not exceed the relative displacement and velocity constraints.

Figure 9 shows the damper force corresponding to the relative velocity. All damper forces are constrained to the bounded region. It can be found that maximum damping is seldom used and only at low relative velocities. The possible reason is that the latching happens at the instant when the device velocity becomes zero or very small. At zero or small relative velocity instants, the hybrid MPC has the same function with latching control, i.e., locking the PTO for certain time intervals to achieve the optimum phase. Furthermore, it can be found that most of the damper force lies on the constraint boundary. So, the boundary coefficients of the damper force significantly affect the hybrid MPC, which will be discussed in the next section.

4.2. Control Force Constraint

In Section 2.4, the method of obtaining variable damping was discussed. It must be noted that the opening aperture of the electronic valves is not unlimited due to hardware size limitations, which also leads to the damper coefficient not being infinitely small, i.e., the parameters β₃ and β₄ cannot be zero. To investigate the effect of parameters β₃ and β₄ on the performance of hybrid MPC, we assumed that β₃ is equal to β₄ and the wave excitation force prediction is exact. The average power $\overline{P}$ over the horizon length N_s in terms of the number of samples can be calculated as follows:

$$\overline{P}=\frac{1}{N_s}{\displaystyle \sum _{t=1}^{N_s}{c}_{pto}{\left[{\dot{x}}_1(t)-{\dot{x}}_2(t)\right]}^2}$$

(39)

The average generated power and power improvement percentage relative to uncontrolled power generation under different parameters β₃ and β₄ are listed in

Table 5. Average power comparison under different parameters β₃ and β₄.

β3, β4 (104)	Average Power (kW)	Power Improvement Percentage (%)
0	82.214	7.4454
0.01	82.206	7.4349
0.1	82.139	7.3474
1	81.469	6.4718
10	75.406	−1.4520
50	55.483	−27.4893
100	41.453	−45.8251

. The simulated result is under the wave condition of T_p = 8 s and H_s = 2.5 m, and the prediction horizon N_p is 32 (8 s). The average generated power of uncontrolled WEC is 76.808 kW.

It can be seen that the average generated power of hybrid MPC with zero β₃ and β₄ is maximum. The power generation gradually decreases as the coefficients β₃ and β₄ increase, and there is a critical value β_c (10⁴ < β_c < 10⁵); the power generation will be less than the case of no control when the boundary coefficients are greater than this value. The variable damper is placed in parallel with the PTO so that the captured energy by WEC is converted into PTO energy and hydraulic energy, with the former being subsequently converted into available electricity and the latter being converted into thermal energy dissipation. The energy converted via internal damping can be improved by adjusting the damper in a small interval of damping coefficient when less energy is dissipated by the variable damper, and most energy is converted into electrical energy. However, a too-large damping coefficient will result in too much energy being dissipated by the damper, reducing the improved performance of variable damping. In reality, a minor boundary coefficient of the electrohydraulic damper is challenging to obtain due to the limitations of the openings, and so a larger boundary coefficient (β₃ = β₄ = 1 × 10⁴ Ns/m) was applied in subsequent studies.

4.3. Prediction Horizon and Wave Condition

The power ratio heatmaps of different prediction horizons under various wave conditions are shown in Figure A1. Each wave state has the same random wave phase. The simulation results are the ratios of generation power of the hybrid MPC system to uncontrolled output power in each wave state. The H_s and T_p of irregular waves range from 1~5 m to 6~12 s, respectively, which cover most of the possible sea conditions.

As can be seen from Figure A1, the hybrid MPC has a negative effect on power generation in a large portion of the wave conditions when the prediction horizon is short N_p ≤ 16 (4 s). With the increase in prediction horizon N_p ≥ 32 (8 s), the power improvement percentage tends to be steady, and the effect of increasing the prediction horizon on power stops being apparent. This indicates that the proposed hybrid MPC has a high requirement for the prediction capability of wave excitation force. However, the simulations here assume that the future wave excitation force is precisely known, which is different from the actual situation. The prediction of the wave excitation force is subject to some error as the prediction horizon increases, which also affects the performance of the hybrid MPC. The effects of the prediction error of wave excitation force will be discussed in the following section.

In the simulations, the maximum power improvement is 22.5% relative to the uncontrolled system (N_p = 48, T_p = 11 s, H_s = 1 m). It should be noted that the simulations all assume a large boundary coefficient (β₃ = β₄ = 1 × 104 Ns/m). If the smaller coefficient can be selected, the power generation by the hybrid MPC system will improve.

4.4. Prediction Error

As mentioned in the above sections, the discussions assumed that future wave excitation force is precisely known but needs to be predicted in reality. Therefore, it is essential to discuss the impact of prediction error on the performance of the hybrid MPC. Two 60th-order (15 s) auto-regressive (AR) models were selected to predict the wave excitation force of the float and reaction section, respectively. The accuracy of the predicted values was described using a goodness of fit (GoF) metric, which is defined as follows.

$$\mathrm{GoF}=\left(1-\frac{\sqrt{{\displaystyle {\sum }_{k=1}^{N_p}{(f_k-{\widehat{f}}_k)}^2}}}{\sqrt{{\displaystyle {\sum }_{k=1}^{N_p}{f}_k{}^2}}}\right)\times 100$$

(40)

where $f_k$ is the real value and ${\widehat{f}}_k$ is the predicted value.

The time history curves of real and predicted wave excitation for the float and reaction section are shown in Figure 10, and their GoF values are summarized in

Table 6. GoF of the AR model and hybrid MPC performance under prediction error.

Np	GoF1 (%)	GoF2 (%)	Power1 (kW)	Power2 (kW)	Power Improvement Percentage (%)
8	99.99	99.99	75.73	77.02	1.70%
16	99.56	99.79	79.48	79.56	1.00%
24	94.31	96.67	81.31	81.06	−0.31%
32	81.62	91.93	81.47	81.44	−0.04%
40	67.28	78.06	81.34	80.57	−0.95%
48	39.58	60.09	81.28	77.50	−4.65%

(H_s = 2.5 m, T_p = 8 s). The GoF₁ and GoF₂ are the prediction accuracy for the float and reaction section, respectively. It can be seen that the GoF is large within a short prediction horizon, and the prediction error gradually increases as the prediction horizon increases. The AR model has sufficient capacity to predict future wave excitation force when the prediction horizon is less than 6 s, with over 94% GoF for the float and reaction section. However, the prediction accuracy decays rapidly when the prediction horizon exceeds 8 s. Another noteworthy phenomenon is that the prediction accuracy for the reaction section is greater than that of the float for the same prediction horizon. The possible reason is that the motion response of the reaction plate is smooth and tends to be more of a stationary time series.

In Table 6, Power₁ and Power₂ are the generation power under the exact prediction and the generation power under the AR prediction, respectively. It can be seen that the power generation under the inaccurate prediction is better than that under exact prediction within a short horizon (N_p ≤ 16). This indicates that the solution of the hybrid MPC is suboptimal, and a more accurate solution method needs to be investigated. As the prediction error increases, the degree of power degradation for non-accurate prediction tends to be extended. In summary, the determination of the prediction horizon is essential, which has a systemic impact on MPC.

5. Conclusions

In this paper, a variable damper was proposed for controlling the dynamics of a two-body WEC with mechanically driven PTO, which introduces logic constraints in the system. A hybrid MPC formulation involving logic variables and continuous variables was constructed to maximize power extraction with the consideration of model limitations. The influence and function of parameters in the proposed method on the two-body WEC have been numerically studied.

The damper valve opening aperture, which imposes a limit on the control force, has an impact on the performance of the proposed hybrid MPC. The larger the maximum opening aperture that can be achieved by the damper valve, the better the performance of the hybrid MPC that could be obtained. The proposed hybrid MPC has a high requirement for the prediction capability of wave excitation force and has the best performance under long T_p and small H_s wave conditions. The prediction horizon has a systemic impact on hybrid MPC performance. A long prediction horizon can improve generation power compared to no-control devices, although a large prediction error accompanies this.

Moreover, despite the realistic feasibility of the proposed controller, there are still some problems that limit the use of hybrid MPC, such as computational efficiency and excitation force prediction accuracy. These can be the subject of future work.