Passive Model Predictive Control on a Two-Body Self-Referenced Point Absorber Wave Energy Converter

: Wave energy is nowadays one of the most promising renewable energy sources; however, wave energy technology has not reached the fully-commercial stage, yet. One key aspect to achieve this goal is to identify an effective control strategy for each selected Wave Energy Converter (WEC), in order to extract the maximum energy from the waves, while respecting the physical constraints of the device. Model Predictive Control (MPC) can inherently satisfy these requirements. Generally, MPC is formulated as a quadratic programming problem with linear constraints (e.g., on position, speed and Power Take-Off (PTO) force). Since, in the most general case, this control technique requires bidirectional power ﬂow between the PTO system and the grid, it has similar characteristics as reactive control. This means that, under some operating conditions, the energy losses may be equivalent, or even larger, than the energy yielded. As many WECs are designed to only allow unidirectional power ﬂow, it is necessary to set nonlinear constraints. This makes the optimization problem signiﬁcantly more expensive in terms of computational time. This work proposes two MPC control strategies applied to a two-body point absorber that address this issue from two different perspectives: (a) adapting the MPC formulation to passive loading strategy; and (b) adapting linear constraints in the MPC in order to only allow an unidirectional power ﬂow. The results show that the two alternative proposals have similar performance in terms of computational time compared to the regular MPC and obtain considerably more power than the linear passive control, thus proving to be a good option for unidirectional PTO systems.


Introduction
Ocean energy and particularly wave energy is nowadays one of the most promising, but still underutilized, renewable energy sources that can be commercially exploited to provide power to the electrical grid [1]. One key aspect to achieve this goal is an efficient control strategy which absorbs maximum energy from the waves, while respecting the physical constraints of each WEC.
Linear control strategies are still the most used for WECs. They have been deeply studied in the WEC literature and provide a reliable, robust and efficient control in most cases. Among them, in reactive control the force applied by the PTO system is the sum of multiple components, proportional to the position, velocity, and acceleration of the WEC, respectively [2]. This strategy typically requires the PTO to be oversized and work both as generator and motor, which can have a negative impact on its efficiency [3]. In reactive control, the optimal velocity profile is a function of the future wave excitation force; thus, this velocity cannot be determined optimally in practice. For this reason, suboptimal control strategies based on causal functions are used [4]. Moreover, reactive control presents drawbacks from the energy perspective, requiring large bidirectional power flows between the WEC and the PTO [5]. In addition, displacements and velocities may be too high and would have to be considered in the mechanical design and managed in the WEC operation.
Passive control is another strategy where the force applied by the PTO is just proportional to the WEC velocity, which allows a unidirectional power exchange, from the WEC to the PTO. Compared to reactive control, passive control presents the following advantages: (a) unidirectional instead of bidirectional power flow means moderate losses; and (b) the power peaks and the amplitude of the oscillations are lower, so the installation becomes cheaper and less critical. The main passive control drawback is the limited power capture compared to reactive, especially for wave profiles far from the natural response of the WEC.
The Model Predictive Control (MPC) is being used with certain success on Wave Energy Converters, because its formulation is optimal in terms of energy extraction and is able to inherently manage physical constraints [6,7]. This technique requires a WEC model and an excitation force prediction within a predefined time interval in order to predict the next control action on the PTO force [8]. This action is obtained optimizing an objective function over the prediction horizon. Generally, MPC on WEC applications is formulated as a quadratic programming (QP) problem with linear constraints, where the objective function is the extracted energy and the constraints are the position, speed and PTO force [9,10]. Since in the most general case this control technique requires bidirectional power flow between the PTO system and the grid, it has similar disadvantages as reactive control: (a) PTO working as generator and motor; (b) increase of PTO losses; and (c) reduction of efficiency.
Few MPC applications to WECs have addressed the problem of restricting the power flow direction from the WEC to the PTO, thus allowing its application to a broader range of existing WECs. The reason for this gap is that allowing this unidirectional energy flow using QP in MPC requires setting nonlinear constraints. This turns the QP optimization with linear constraints into a Nonlinear Programming (NLP) problem with nonlinear constraints, making the computational effort significantly higher [6]. Among the few contributions addressing the problem in the available literature, O'Sullivan and Lightbody [11] restricted the power flow direction by means of a nonlinear constraint function that depends on the speed, flux, resistance and current of an electrical PTO. To allow just unidirectional power flow, Paparella and Ringwood [12,13] proposed a MPC-like strategy using pseudo-spectral control with truncated Fourier series as orthogonal base to handle the nonlinear power constraints in a more efficient way. Both proposals have, however, the same drawback: nonlinear constraints imply higher processing time to solve the optimization problem. Tom and Yeung [14] studied a permanent magnet generator as PTO with a passive non-linear MPC control strategy. They defined the PTO force as a linear function of the velocity (F pto = b g v(t)) and integrated the dynamic of the WEC inside the objective function from a numerical perspective using the Berkeley Library for Optimization Modeling software (BLOM) on a Matlab/Simulink (version 2012, Berkeley, CA, USA) model. This work proposes two independent MPC control strategies applied to a two-body self-reference point absorber that address this issue from two different perspectives: (a) adapting the MPC formulation to achieve passive loading strategy similar to the work by Tom and Yeung [14], but integrating analytically the dynamic of the WEC inside the objective function; and (b) re-formulating the linear constraints in the MPC in order to only allow an unidirectional energy flow. The results show that the proposed control strategies have similar performance compared to the regular MPC in terms of computational time, and they extract considerably more power than the classical passive control.
This rest of the paper is divided into six sections. Section 2 presents the model of the two-body self-reference point absorber. Section 3 shows the traditional MPC approach considering linear and nonlinear constraints. The two MPC proposals based on passivity (to prevent bidirectional energy flow) are formulated in Sections 4 and 5, respectively. Section 6 presents the analysis of the results. Section 7 draws the conclusions.

Model of the Two-Body Self-Reference Point Absorber WEC
The WEC model used in this work, as shown in Figure 1, was developed within the IMAGINE project [15] funded by the European Union's Horizon 2020 programme. This WEC is a two-body self-referenced point absorber formed by an external floater ring that moves with respect to an internal cylindrical spar attached to the seabed by means of a mooring. This relative movement is exploited by the PTO to extract power from the WEC. The PTO system is a direct-drive device in which a recirculating ball screw is coupled with a three-phase Permanent Magnet Generator (PMG), as shown in Figure 2 [16] and described in [17]. The linear motion of the ball screw is converted into the nut rotary motion. Being the ball nut connected to the rotor of the permanent magnet generator, its mechanical rotational energy is directly converted in electricity.  Linear wave theory was used to model the wave profile. According to the Cumming's formulation, the dynamic Equations (1) and (2) describe the two-body wave energy converter motion in heave mode [18]: where ρ is the water density, g is the gravity acceleration, k m is the mooring stiffness, b f r is the PTO friction coefficient and F pto (t) is the PTO force. For body i, k hi = ρgS i is the hydrostatic coefficient, m i is the mass, S i is the body transversal area, b vi is the viscosity coefficient, A ∞ ii is the added mass at infinite frequency, F ri (t) is the radiation force, F ei (t) is the excitation force and z i (t),ż i (t) andz i (t) are the position, speed, and acceleration. In the paragraph above, the subscript i is 1 when referred to Body 1 (floater) and 2 when referred to Body 2 (spar). In addition, the PTO mass is considered negligible compared to m 1 and m 2 .

Space State Model in Continuous Time
Both radiation forces F r1 (t) and F r2 (t) above can be modeled by means of state space equations based on the work of Falnes and Yu [19]. The space state models used for the floater and the spar are the following third-and second-order systems, respectively: where where The states of the radiation forces F r1 (t) and F r2 (t) are U 1 (t) = [ u 11 (t) u 21 (t) u 31 (t) ] T and U 2 (t) = [ u 12 (t) u 22 (t) ] T , respectively. The dynamics described by the system of differential equations from (1) to (10) can be modeled by means of the following state space equations [18,19]: where and is the position of the PTO andż 1 (t) −ż 2 (t) is its velocity. Numerical values of the WEC parameters present from Equation (3) to (18) can be found in Appendix A.

Space State Model in Discrete Time
To design an MPC strategy, the discretization of the space state model is necessary. In this section, two space state models in discrete time are presented. The first one uses the discretized PTO force, F pto (k) (where k is the discretized time variable), as control signal and the second uses the variation of the discretized PTO force between successive time instants, ∆F pto (k + 1) = F pto (k + 1) − F pto (k) including the PTO force, F pto (k), as a system state. Both approaches are equivalent and are used as needed by the proposals. The first approach makes it easy to rewrite the control variable F pto (k) by b pto (ż 1 (t) −ż 2 (t)) in order to formulate a first passive MPC strategy as a function of the b pto parameter. The second MPC approach allows constraining the variation of the PTO force between time steps and at the same time assign the linear constraints on the state variables speed and PTO force.

First Approach: PTO Force F pto (k) as Control Signal
Using the discrete solution of the state space in Equation (11) described by Chen [20] and considering the variables constant between samples, in the zero-order hold approximation, the following discrete state space equation is obtained: where is the discretized PTO force, F e1 (k) and F e2 (k) are the discretized excitation forces on both floater and spar and X(k) is the discretized state vector:  (19) including the excitation forces and the PTO force as states. In addition, incremental values of the excitation forces ∆F e1 (k + 1), ∆F e2 (k + 1) and PTO force ∆F pto (k + 1) are the new inputs (Equations (22)-(24)).
The state space equations in terms of the new augmented state vector are: where and 0 x,y is a null elements matrix with x rows and y columns. Setting the incremental values of the PTO force ∆F pto (k + 1) as a control signal allows avoiding abrupt changes in the PTO force between consecutive control actions using an additional linear constraint in the MPC optimization problem. This new constraint is explained in Section 3.3.

Traditional Model Predictive Control
MPC gets an optimal control for optimizing an objective function (minimization or maximization) over a finite future time horizon. To implement an MPC strategy successfully, the following are necessary [22]: (a) a mathematical model of the system to predict the output on the future horizon; (b) an objective function to obtain an optimal control series over the future horizon using an optimization procedure with constraints; and (c) a receding strategy where, at each instant, the first sample of the optimal control series calculated by the optimization is applied.

Prediction Equations of the Traditional MPC
To get the prediction equations, the methodology shown by Rossiter [23] is used. This methodology is based on the discrete state space in Equation (25) obtained by assuming the variation of the PTO force ∆F pto (k) between two consecutive time instants as the control signal (second approach in Section 2.2): and n h is the sample number of the prediction horizon.

Objective Function of the Traditional MPC Control
The objective function J(k) to be minimized at each time instant, k, is the opposite of the average power absorbed by the PTO over the prediction horizon.
The objective function above can be rewritten in terms of the future output V − → (k) of Equation (27).
where Q 1 is the following square matrix with order 3n h .
Finally, substituting Equation (27) into Equation (35), the following MPC objective function of a typical quadratic programming problem is obtained:

Constraints on the Objective Function for the Traditional MPC
To set maximum and minimum values to the force, velocity, position and force increment of the PTO, linear constraints on the MPC objective function can be added to the optimization problem [7]. Restrictions on the PTO force increments can be simply formulated as: On the other hand, constraints on the position, speed and PTO force can be set in the prediction vector V(k) (Equation (29)) and expressed by means of Equation (27) as a linear matrix relationship as a function of the control variable ∆F pto (k): Regarding the restrictions defined in Equations (39) and (40), the optimization problem of Equation (38) becomes a quadratic programming problem with linear constraints.

Non-Linear Constraints on the Objective Function for the Traditional MPC to Get Passivity
A direct way to get unidirectional power flow in the MPC optimization problem of the Equation (38) is constraining the PTO power at each instant by means of the product between PTO force and velocity [11]. However, as mentioned in Section 1, this approach is expensive in terms of computational time, hence it provides a solution to the passivity problem that is not viable for real-time applications, thus limiting its practical usability. In the following, the PTO power constraint is enforced at each time instant, k, by means of the product between PTO force and velocity: The MPC objective function of Equation (38) and the non-linear constraint of Equation (41) define a nonlinear programming problem (NLP) with nonlinear constraints.

Passive MPC Proposal 1: Adapting MPC Formulation to Passive Loading Strategy
One of the passive control strategies is the resistive load control where the PTO force just depends on the PTO velocity in the following linear way: To constraint the PTO force to take values such that the control is passive, Equation (42) is discretized and substituted in the discrete model of Equation (19), which was obtained by using the PTO force F pto (k) as the control signal (first approach in Section 2.2): The last term in Equation (43) can be expressed using the state space vector X(k). From Equation (21):ż where Finally, the discrete space state model can be formulated as: where Notice that matrix A d2 is function of b pto . Equations (46) and (47) model a two-body point absorber with passive control in discrete time, where the PTO damping is defined by b pto .

Prediction Equations of the Passive MPC Proposal 1
The main idea of this passive MPC proposal is to determine the PTO damping at each time instant, k, for the WEC modeled by Equations (46) and (47). This means that: To do so, an optimization problem is formulated using the discrete system model of Equations (46) and (47) and prediction of the excitation forces on a future horizon. Based on the discrete state space in Equations (46) and (47) and using the methodology shown by Rossiter [23], the following prediction equations for the passive MPC formulation can be obtained: where . . .

C d A d2
Notice that the matrices P, H e1 and H e2 in Equations (50)-(53) depend on A d2 , which is function of b pto (k) (Equations (48)-(49)). It means that these matrices change in each time instant, k, as b pto (k) changes.

Objective Function of the Passive MPC Proposal 1
In passive control, the opposite of the average power absorbed at each time instant, k, by the PTO over the prediction horizon can be formulated as: Equation (54) can be rewritten in term of the future output Y − → (k) from Equation (50).
where Q 2 is the following square matrix with order 2n h .

Constraints on the Objective Function of the Passive MPC Proposal 1
Although there is no constraint in the optimization of the objective function given in Equation (58), the positive sign of b pto (k) at each time instant k should be verified.

Passive MPC Proposal 2: Adapting Linear Constraints in the MPC to Get Passivity
This passive MPC proposal uses the same objective function as the traditional MPC (Equation (38)). Unlike the first proposal in Section 4, which determines PTO damping b pto (k) at each time instant, this MPC optimization problem leads to estimate a PTO force series over the prediction horizon at each time instant. Passive control is imposed by properly managing the linear constraints defined in Equations (39) and (40). The idea is to define speed and PTO force constraints at each time instant over the prediction horizon in such way that the product between speed and PTO force is positive.

Prediction Equations of the Passive MPC Proposal 2
This proposal uses the same prediction Equation (27) as the traditional MPC.

Objective Function of the Passive MPC Proposal 2
This proposal uses the same objective function Equation (38) as the traditional MPC.

Constraints on the Objective Function of the Passive MPC Proposal 2
Linear constraints in the traditional MPC (Equation (40)) allow estimating PTO force increments ∆F pto − −− → (k) over the prediction horizon, limiting the variation range of the position (x min pto , x max pto ) , velocity (v min pto , v max pto ) and force of the PTO (F min pto , F max pto ) by setting maximum and minimum values in the vector V − → (k) (Equation (29)). Generally, these limits are constants over the prediction horizon and for all time instants k, hence, the elements of the vectors V − → (k) min and V − → (k) max defined in Equation (40) are: To set a passive MPC strategy that handles the linear constraints, maximum and minimum values of Equation (59) must be considered as variants over the prediction horizon and for every time instant k: The idea is to find values to the elements of the matrices in Equations (60)  Linear constraints of the MPC are assigned at each time instant (k + i) over H: This way of constraining velocity and force of the PTO will guarantee unidirectional power flow from the WEC to the PTO when MPC is applied.

Numerical Analysis
The WEC of Figure 1 is a the two-body self-referenced point absorber formed by an external floater ring with a 19.2 m diameter and 4.8 m height that moves with respect to an internal cylindrical spar with a 5.7 m diameter and 36.5 m height attached to the seabed by means of a mooring. This WEC is inspired by the Ocean Power Technology (OPT) PowerBuoy TM . Similar two-body WECs are researched in [24] as RM3 and in the NumWEC project [25] as F-2HB. The two-body WEC model used in this work is described in more detail in Appendix A.
The time-domain simulation model of the two-body WEC shown in Figure 1 was implemented in Matlab/Simulink to test the proposed control strategies. The model is formed by three main subsystems: the first subsystem generates the excitation force, the second one models the WEC using a state space equations and the last one implements the control strategy and sets the PTO force.
The Wamit program was used to calculate the hydrodynamic parameters as functions of the frequency. The irregular wave profiles were generated randomly using the Bretschneider spectrum. For each peak period considered, a 15 min simulation was executed. The time step used in all MPC strategies is 0.05 s and the prediction horizon is 9 s. In addition, to focus on the control strategy [11], it is assumed that the future excitation wave forces are perfectly known (exact prediction of the excitation forces) over the prediction horizon.
The results applying the passive MPC Proposal 1 of Section 4 (P-MPC-1) and the passive MPC Proposal 2 of Section 5 (P-MPC-2) were compared with those obtained through the traditional MPC without constraints of Section 3.2 (MPC), passive MPC with nonlinear constraint (P-MPC-NLC) of Section 3.4 (Equation (41)) and passive (i.e., resistive) loading control (RL).
In addition, the MPC strategy with relaxed linear constraints of Section 3.3 (Equations (39) and (40)) was included to consider the computational effort when linear restrictions are added (MPC-LC). In this case, the relaxed constraints are: maximum excursion of 5 m, maximum speed of 5 m/s and maximum PTO force of 10 MN. It is important to highlight that these relaxed constraints are not achieved during the simulations and are only considered to estimate the processing time when linear constraints are added to the MPC.
In addition, an optimal passive strategy (ORL) with PTO constant damping was included as reference control where the best PTO damping in terms of power extracted is obtained by means of an optimization process on the same wave profile. These strategies are explained briefly in Appendix B.
Taking into account that P-MPC-1 had good performances compared to RL and ORL strategies, it was used in Step 3 of the algorithm in Section 5.3 as starting point of P-MPC-2. Further, the constants α and β in Step 4 of the same algorithm were set to 3 and 2.25, respectively, and the excursion constraints were set to x max pto = 2.5 m and x min pto = −2.5 m. The selection of α and β affects the processing time of P-MPC-2 and it was performed to set constraints for the speed and the PTO force that are non achieved, but close to the maximum obtained.
To conveniently compare the control strategies, the maximum number of iterations in the optimization problem associated to the control P-MPC-2 and P-MPC-NLC was adjusted to obtain processing times in the same range as the control MPC-LC (less than 1 s). In this sense, the number of iterations in the control P-MPC-2 and P-MPC-NLC were 5 and 50, respectively. Besides, for the other control strategies, the maximum number of iterations used was the value set by default (250) in the Matlab function of the optimization problem.

Relevant Average and Maximum Results
This section presents the relevant results considering irregular waves, varying the peak period from 5 to 16 s and considering significant wave heights of 1.25, 3.25 and 5.25 m. The seed to create each wave profile was changed. This means that the wave profiles are not the scaled version of each other. Figure 3 shows the average PTO power for all considered wave profiles and for the control strategies RL, ORL, MPC, P-MPC-1, P-MPC-2, P-MPC-NLC and MPC-LC described above. It can be noted that the traditional MPC strategies (MPC and MPC-LC) allow higher average power absorption due to their reactive nature. This means that the energy is circulating in a bidirectional way between the PTO and the WEC.
On the other hand, the control strategy P-MPC-2 presents better performances than P-MPC-1 in terms of extracted power. The reason is that P-MPC-1 is a simple and faster strategy that keeps the same PTO damping over the entire prediction horizon, while P-MPC-2 has to estimate the PTO force for each time instant over the prediction horizon being more complex and slower.
Taking into account the differences in processing time and the power extracted between both proposed control strategies P-MPC-1 and P-MPC-2, it is convenient to divide the analysis in two groups. The first group is composed by faster control strategies with lower extracted power (P-MPC-1, RL and ORL), and the second one is composed by slower control strategies with higher extracted power (P-MPC-2 and P-MPC-NLC). A third group formed by the reactive control strategies with the greatest extracted power (MPC and MPC-LC) serves as reference for the first and second groups. The MPC strategy without constraints is a convex quadratic programming problem whose solution is fast in terms of processing time. For this reason, the first group is compared to the MPC strategy. Besides, the MPC-LC strategy considers the computational effort when linear restrictions are added. These linear constraints are relaxed, never achieved and added to estimate the processing time in the quadratic programming problem with linear constraints (non-convex programming problem). For this reason, the second group, formed by control strategies, for which passivity is defined in their constraints, is compared to the MPC-LC strategy.
Based on Figure 3, it can be concluded that the proposal P-MPC-1 extracts 5% more average power than the RL strategy and practically the same average power as the ORL reference strategy. On the other hand, the proposal P-MPC-2 allows obtaining 4% more power than the strategy P-MPC-NLC. Concerning the average processing time to solve the associated optimization problem, Figure 4 shows that the most time-consuming control strategy is P-MPC-NLC (0.973 s), followed in descending order by MPC-LC (0.844 s), P-MPC-2 (0.661 s), P-MPC-1 (0.029 s) and MPC (0.013 s). It can be highlighted that the processing time of P-MPC-1 is low and in the same order as MPC. On the other hand, the processing time of the control strategy P-MPC-2 is 32% lower than P-MPC-NLC and 22% lower than MPC-LC. According to the results reported in Figures 3 and 4, both proposals are competitive in terms of absorbed power and processing time. Figure 5 shows that the maximum speed takes high values for MPC and P-MPC-NLC, low values for the control RL and ORL and intermediate values for both proposals P-MPC-1 and P-MPC-2.       Figure 9 presents the values of the PTO damping as a function of time. It can be noted that the values in both P-MPC-1 and P-MPC-2 vary over time while the PTO damping in the RL and ORL strategies is constant. In addition, the PTO damping in P-MPC-1 presents less abrupt changes and lower values compared to P-MPC-2. The reason for this is that P-MPC-1 is based on keeping the PTO damping constant over the prediction horizon, while P-MPC-2 estimates a different PTO damping for each time instant over the prediction horizon.   To vary the processing time, the number of iterations was increased conveniently on both strategies. The wave profile used in this analysis is the sea state with peak period Tp = 12 s and significant wave height Hs = 3.25 m.

Sensitivity Analysis
In Figure 10, it can be noted that a higher average power is obtained applying P-MPC-2, when the average time spent in the optimization problem is less than 1.15 s., while P-MPC-NLC extracted more power for a higher processing times (greater than 1.15 s).
Since real-time control requires fast response from the processors in order to get an estimation of the control action, the control strategy P-MPC-NLC is not feasible from a technical point of view. On the other hand, as the proposal P-MPC-2 consumed less processing time extracting more average power, it could be suitable for real time implementation.

Conclusions
This work proposes two passive MPC control strategies applied to a two-body selfreference point absorber. The first proposal adapts the MPC formulation to a passive loading strategy (P-MPC-1) and the second one adapts the linear constraints in the MPC in order to allow only unidirectional power flow between the wave energy converter and the power take-off (P-MPC-2).
The results show that the first proposal has a processing time in the same order (0.029 s) as the MPC without constraints (0.013 s), extracting 5% more power than the resistive loading (RL) and 0.1% more power than the passive reference control strategy (ORL), achieving intermediate values for PTO speed and low values for PTO force.
On the other hand, the second proposal obtains 4% more average power than the passive MPC strategy with non-linear constraints, consuming 32% less processing time than the passive MPC strategy with non-linear constraints and 22% less processing time than MPC with linear constraints. In addition, this proposal presents intermediate values to the maximum speed and the PTO force.
Due to its low processing time, the first proposal P-MPC-1 could be used to improve the performance of the traditional passive control (RL) including constraint. In this sense, in this article, the first proposal is used to initialize the second proposal. On the other hand, the second proposal could be used to implement a technically feasible passive MPC.
As an MPC strategy, our proposals are very dependent on the model used in their formulation as well as on the prediction of the excitation force. A non-representative model of the WEC would lead to an inefficient MPC strategy.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
The WEC of Figure 1  Both radiation forces Fr 1 (t) and Fr 2 (t) were modeled by means of state space equations using third-and second-degree systems: The dynamics of the system was modeled by means of the following state space equations:Ẋ (t) = AX(t) + B e1 F e1 (t) + B e2 F e2 (t) + B pto F pto (t)