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

: In this paper, a variable damper is proposed to regulate the efﬁciency 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 signiﬁcant wave height and long wave period, improving the power generation of the uncontrolled system up to 22.5%. And, the prediction error has a signiﬁcant effect on hybrid MPC performance, especially for long prediction horizon.


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 twobody 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.

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]: .. ..
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. Figure 1. Overview of the two-body wave energy converter system.

Hydrodynamic Force
The regular wave excitation force can be represented as follows: where  denotes the real part of the formula; H is the wave height; and ( ) exc F  is the frequency-dependent complex amplitude of the wave excitation force. The dependence of ( ) exc F  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: where Aj is the wave amplitude at angular frequency j  , which can be obtained from the wave spectrum; j  is the randomized phase angle at angular frequency j  , and N is the number of frequency bands selected to discretize the wave spectrum. The Pierson-Moskowi (PM) spectrum was chosen to describe incident waves in this paper.
where f is the frequency of the PM spectrum; Hs is the significant wave height; fp is the peak frequency; and dωj is the frequency interval. The radiation force can be represented as follows [23]:

Hydrodynamic Force
The regular wave excitation force can be represented as follows: where denotes the real part of the formula; H is the wave height; and F exc (ω) is the frequency-dependent complex amplitude of the wave excitation force. The dependence of F exc (ω) 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: where A j is the wave amplitude at angular frequency ω j , which can be obtained from the wave spectrum; φ j is the randomized phase angle at angular frequency ω 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.
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]: x j (τ)dτ (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: 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 . 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 × 2. where ij is the added mass at infinite frequency and Krij 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: where Ar is the system matrix, the dimension of which depends on the accuracy of approximation representation; Br is the input matrix; Cr is the output matrix; Dr is the feedthrough matrix, which is assumed to be zero to maintain the causality of the system and will be omi ed in the following equations; and q and 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 2  . The net buoyancy restoring force can be represented as follows: where ρ is the seawater density; g is the gravity acceleration; S is the water plane area of the body; and KH is the hydrostatic stiffness coefficient.

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: where ω is the ballscrew angular speed, that is, the input angular speed of the gearbox; 1 2 x x    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: where Mi is the input moment of the gearbox, and the transmission ratio of the gearbox is denoted as ng. The angular speed and input moment of the generator are calculated as follows: The net buoyancy restoring force can be represented as follows: 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.

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: where ω is the ballscrew angular speed, that is, the input angular speed of the gearbox; .
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: 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: 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.
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.
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].
The PTO mechanism can be represented as a linear mass-spring-damper system, where the PTO force is given by the following equation.
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.

Control Force
In the proposed controller, the control force changes with the variable damping, and the damper force can be represented as follows: 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 . Compared to the classical passive hydraulic damper, the electrohydraulic d comprises electronic valves instead of passive valves. The required damping coe can be achieved by adjusting the opening aperture of electronic valves appropriatel trolled by the current i of the servo mechanisms. The energy consumed in obtain appropriate damping coefficient is relatively small or negligible. A typical force-c map for an electrohydraulic damper is shown in Figure 4. It should be noted that th itive direction of damper force Fd is opposite to that of the control force u F , i.e., Fd  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.
The required control current for driving the variable damper system is based on calculated control force Fu and relative velocity 1 2 x x    at each time, which can be described as a non-linear inverse function of the damper force, as follows: 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 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.
The required control current for driving the variable damper system is based on calculated control force F u and relative velocity .
x 2 at each time, which can be described as a non-linear inverse function of the damper force, as follows: 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.

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.
where F c rad1 and F c rad2 are the radiation forces (convolution integral term) loading in the float and reaction section, respectively. q ij , i, j ∈ {1, 2} is the system state of the state-space realization for the convolution term of radiation forces F radij .
Substituting Equations (6), (8), (15) and (21) into Equation (1), the dynamic equation can be restated as follows: where, It can be assumed that the state vector is chosen as x 2 q T T ; then, the continuous state-space model can be obtained by substituting Equation (20) into Equation (22). 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. where,

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 .
where the connective "↔" means "if and only if". This mixed logic continuous relation can be translated into linear constraints involving the logical variable δ v . .
where v l and v u are lower and upper bounds to the function .
The constraints for damper force can be written in a more detailed form, as follows [28]: where z 1 , z 2 , and z 3 ∈ R are continuous auxiliary variables, which are imposed constraints.
A generic if-then-else construct of the form is as follows: where δ ∈ {0, 1}, z ∈ R, x ∈ R n , u ∈ R m , a 1 , a 2 , b 1 , b 2 , f 1 , and f 2 are constants with approximate dimensions, which can be translated into integer linear inequalities of the form.
where f i and F i are lower and upper bounds to a i x + b i u + c i , i ∈ {1, 2}. 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. x where x min and x max are the maximum and minimum relative displacement that the system can achieve. .
x min and .
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.
In this paper, p ∈ R 12 is the state vector; F d ∈ R is the control input;δ v ∈ {0, 1} is the logic variable; z ∈ R 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 2 , E 3 , E 1 , E 4 , and E 5 are matrices of suitable dimensions.

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: For convenience, the objective function can be written as a cost function (and thus is a minimization problem) of the form.
where N p is the optimization horizon length.
Based on the hybrid logic constraints, the hybrid MPC can be formulated as follows:

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 (the parameters are from Reference Model 3 [31]). The seawater density ρ takes the value of 1025 kg/m 3 and the gravity acceleration takes the value of 9.81 m/s 2 . 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. The boundary coefficients of damper force are shown in Table 3, and the relative displacement and velocity constraints are summarized in Table 4.   The same PTO mechanism was used for all simulations. The parameters of the PTO system can be found in Table 2. The boundary coefficients of damper force are shown in Table 3, and the relative displacement and velocity constraints are summarized in Table 4.   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 Figures 6 and 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 o 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 zer  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 forces are constrained to the bounded region. It can be found that maximum dam seldom used and only at low relative velocities. The possible reason is that the happens at the instant when the device velocity becomes zero or very small. At small relative velocity instants, the hybrid MPC has the same function with latch trol, i.e., locking the PTO for certain time intervals to achieve the optimum phase. more, it can be found that most of the damper force lies on the constraint bound the boundary coefficients of the damper force significantly affect the hybrid MPC will be discussed in the next section.  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.

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 β3 and β4 cannot be zero. To investigate the effect of parameters β3 and β4 on the performance of hybrid MPC, we assumed that β3 is equal to β4 and the wave excitation force prediction is exact. The average power P over the horizon length Ns in terms

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 β 3 and β 4 cannot be zero. To investigate the effect of parameters β 3 and β 4 on the performance of hybrid MPC, we assumed that β 3 is equal to β 4 and the wave excitation force prediction is exact. The average power P over the horizon length N s in terms of the number of samples can be calculated as follows: The average generated power and power improvement percentage relative to uncontrolled power generation under different parameters β 3 and β 4 are listed in Table 5. 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 β 3 and β 4 is maximum. The power generation gradually decreases as the coefficients β 3 and β 4 increase, and there is a critical value β c (10 4 < β c < 10 5 ); 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 (β 3 = β 4 = 1 × 10 4 Ns/m) was applied in subsequent studies.

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 (β 3 = β 4 = 1 × 104 Ns/m). If the smaller coefficient can be selected, the power generation by the hybrid MPC system will improve.

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.
where f k is the real value andf 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 (H s = 2.5 m, T p = 8 s). The GoF 1 and GoF 2 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 1 and Power 2 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.  In Table 6, Power1 and Power2 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 be er than that under exact prediction within a short horizon (Np ≤ 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.

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.

Appendix A
for the prediction capability of wave excitation force and has the best performance under long Tp and small Hs 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.     Figure A1. Average output power ratio of different prediction horizons under irregular wave conditions.