A Study about Performance and Robustness of Model Predictive Controllers in a WEC System

Abstract

This work is located in a growing sector within the field of renewable energies, wave energy converters (WECs). Specifically, it focuses on one of the point absorber waves (PAWs) of the hybrid platform W2POWER. With the aim of maximizing the mechanical power extracted from the waves by these WECs and reducing their mechanical fatigue, the design of five different model predictive controllers (MPCs) with hard and soft constraints has been carried out. As a contribution of this paper, two of the MPCs have been designed with the addition of an embedded integrator. In order to analyze and compare the MPCs with conventional PI type control, an exhaustive study about performance and robustness is realized through the computer simulations carried out, in which uncertainties in the WEC dynamics and JONSWAP spectrum are considered. The results obtained show how the MPCs with embedded integrator improve power production of the WEC system studied in this work.

1. Introduction

Nowadays, the main motivation for the research and development of wave energy converters (WECs) is the advantages offered by waves: a clean and abundant source of energy. As evidence, the authors in [1] compared a global study of net wave power (estimated at about 3 TW) with the electrical power consumed globally in 2008 (equivalent to an average power of 2.3 TW). However, it should be noted that currently, there is not a clear line of development, but a great diversity of systems based on different approaches to extract energy from the waves. In particular, the ocean energy systems collaboration program [2] classifies three kinds of WEC systems: oscillating water columns, overtopping and wave-activated bodies (WABs). This work focuses on a type of WAB system, a point absorber wave (PAW) energy converter from the W2POWER (Wind and Wave Power) platform [3]. These systems are characterized because their extension is significantly smaller than the predominant wavelengths. In addition, PAWs extract the maximum mechanical power from the sea when they are in resonance with the excitation force caused by the waves [4]. In order to favor this situation, it is necessary to enlarge the bandwidth of these devices. For this reason, several control systems are used, and the most common are: passive loading control, reactive loading control and latching control [5,6]. Although, since the last decade, more complex controllers like MPC (model predictive control) are being employed. The interest in implementing MPCs in WEC systems is motivated by the need to increase the productive/economic viability of these systems; due to the fact that these controllers allow them to minimize mechanical fatigue in their structures (limiting the operating ranges) and to focus the control strategy on the maximization of the extracted power directly.

Actually, several authors have developed predictive controllers for WECs [7,8,9,10,11,12,13,14,15,16,17]. These MPCs can be grouped according to the characteristics of the cost function optimized. On the one hand, the authors in [8,10,11,12,13,14,15,16] used a cost function in which the extraction of mechanical power was directly maximized. Whereas, on the other hand, in [9,17], the authors proposed a cost function to maximize power extraction by minimizing the error between the speed of oscillation of the system and a setpoint for it. In addition to the above classification, MPCs can be distinguished according to the mathematical model used for their design. On one side, in [9,10,11,13,14], a reduced model was used for the design, which did not consider the dynamics of the radiation force. Meanwhile, in [7,8,12,13,15,16,17], such dynamics were considered in the design model. Moreover, the previous works did not consider the dynamics of the power take-off system (PTO), neither in the evaluation, nor in the design of the MPCs. In this aspect, the authors in [11,14], although they did not consider the PTO dynamics, optimized the cost function for the increment of the control signal, thus limiting the slew rate of the actuator. Finally, the treatment carried out in the case of non-feasibility when solving the cost functions with restrictions should be highlighted. In this aspect, the works in [9,10,13] considered a more complete approach by adding soft constrains in the case of non-feasibility.

This paper analyses the main approaches of MPCs applied to PAW systems [8,9,10,11,12,13,14,15,16,17]. In addition, a new design is proposed: MPC based on a model with an embedded integrator for controllers that follow a setpoint for the speed of oscillation of the PAW. This approach is recommended in the theory of predictive control in the space of states [18,19], and it is proposed as an alternative method to the one carried out in [9,17] for PAW systems. Moreover, all predictive controllers of this work consider soft and hard constraints. On the other hand, the PTO dynamics is taken into account to validate the MPC controllers, as well as in the design of some of these controllers. After an exhaustive fine-tuning for all MPCs, the main contribution of this work is obtained, and an in-depth study about performance and robustness of all MPCs through simulations is carried out. The results obtained for a sea state defined by the JONSWAP spectrum [15,16] are compared with conventional controllers: I-P control and resistive damping (RD). Furthermore, the mathematical model of the WEC system is obtained using a software simulation based on the boundary element method, such as openWEC [20].

The rest of the article is organized as follows: Section 2 presents the generic mathematical model of a PAW, the standard identification methodology applied in this paper and the treatment applied to the model identified for its later use in the design of the MPCs. Section 3 details the five MPC controllers designed with hard and soft restrictions. Section 4 presents a performance and robustness comparison analysis for the MPCs and conventional (P, PI) controllers. Finally, in Section 5, the main conclusions are indicated.

2. Mathematical Model

Given the importance of mathematical modeling in the design of MPC controllers, this section begins by describing the generic model of a PAW. This is followed by the standard identification process realized in this paper for the study system, one of the WEC systems of the platform W2POWER [3]. Later, the treatment of the model is detailed for its later use in the design of MPC controllers.

2.1. Generic Mathematical Model for Point Absorber Wave

The modeling of the forces affecting a PAW that extracts power from the waves using a single degree of freedom (heave) is widely used in the literature [4,7,8,9,21,22,23,24,25,26,27]; see Figure 1. The main dynamic interactions between the buoy and the waves are collected by:

$$m\ddot{z}=F_e+F_r+F_s+F_u$$

(1)

where m is the mass of the system, z the vertical displacement of the buoy, $F_e$ the excitation force caused by the wave, $F_r$ the radiation force, $F_s$ the hydrostatic restoring force and $F_u$ the control force realized by the PTO system.

The force of radiation is due to the effect produced on the system by the waves it radiates when oscillating. It is modeled by the Cummins equation [28] as (2), where the radiation force $F_r$ is composed of two terms: $F{r}_{m_{\infty }}$ and $F{r}_{K_r}$. The first is a function of the acceleration of the system and the mass of water added $m_{\infty }$; while the second defines a radiation force as a function of the speed of oscillation as an integral of convolution.

$$F_r=-\underset{F{r}_{m_{\infty }}}{\underbrace{m_{\infty }\dot{w}(t)}}-\underset{F{r}_{K_r}}{\underbrace{{\int }_{\infty }^t{h}_r(t-\tau )w(\tau )d\tau }}$$

(2)

The convolution term represents the impulse response that relates the speed of the oscillation system with the force of radiation, where $h_r(t)$ represents the radiation impulse response function (RIRF). To avoid convolution calculations [21,24,26], an approximation is made based on ordinary differential equations (ODE) or as a transfer function in the Laplace domain (3).

$${\displaystyle \frac{F{r}_{K_r}(s)}{W(s)}}=K_r(s)={\displaystyle \frac{a_n{s}^n+a_{n-1}{s}^{n-1}+\dots +a_0}{b_n{s}^n+b_{n-1}{s}^{n-1}+\dots +b_0}}$$

(3)

The hydrostatic restoring force represents the effect of Archimedes and gravity on the buoy (4). By linearizing Equation (4), the hydrostatic force is approximated as (5).

$$F_{res}=-(V_{desp}(z)\rho g-mg)$$

(4)

$$\begin{array}{cc}\hfill F_{res}& =-k_{res}z(t)\hfill \\ \hfill k_{res}& =\rho g{A}_W\hfill \end{array}$$

(5)

where $V_{desp}$ is the volume of water displaced, $\rho$ is the density of seawater, g is the gravity constant, $k_{res}$ is the hydrostatic restoring coefficient and $A_W$ is the area of the buoy on its waterline.

For the excitation force caused by the waves, the components with the highest frequency are negligible. For this reason, authors such as [7,8,11,12,15,16,24,25,26,27] have defined the excitation force as a low pass filter of first/second order at wave height (6). This modeling is used later in the performance study.

$$G_{Fe}(s)={\displaystyle \frac{Fe(s)}{\eta (s)}}={\displaystyle \frac{K_{\tau }}{s^2+2{\zeta }_{\tau }{\omega }_{n\tau }s+{\omega }_{n\tau }^2}}.$$

(6)

On the other hand, in order to model the force on a buoy, the Morison equation [29] can be used, which is habitually employed to estimate the wave loads in the design of offshore structures. By linearizing the Morison equation for a point absorber wave in heave, the excitation force is defined according to Equation (7). This modeling is used later in the robustness study.

$$F_e(t)=m_{\infty }\ddot{\eta }(t)+B_{aprox}\dot{\eta }(t)+k_{res}\eta (t)$$

(7)

where $\eta$ is the height of the incoming wave and $B_{aprox}$ represents the damping of the system, which can be obtained as the stationary gain of the transfer function (3).

The PTO dynamics can be approximated as a linear second order system with a force limitation. The linear relation between the demanded force by the control system, ($F_u$), and the force applied to the buoy, ($F_{pto}$), is given by (8).

$$G_{pto}(s)={\displaystyle \frac{F_{pto}(s)}{F_u(s)}}={\displaystyle \frac{d_a{{\omega }_n}^2}{s^2+2\zeta {\omega }_{n}s+{{\omega }_n}^2}},$$

(8)

2.2. Forces Identification Using Simulation Based on BEM

A common approach to determine the hydrodynamic forces of interaction between the wave and buoy is to use the linear wave theory, which assumes that waves are the sum of incident, radiated and diffracted wave components [21]. These components can be modeled using linear coefficients obtained from the frequency-domain potential flow BEM (boundary element method) solver Nemoh [30]. The BEM solutions are obtained by solving the Laplace equation for the velocity potential, which assumes that the flow is inviscid, incompressible and irrotational. In this work, openWEC software is employed. It is an open-source tool to simulate the hydrodynamic behavior and energy yield from single-body wave energy converters [20]. Two software packages are coupled in openWEC, a frequency domain solver Nemoh and a time domain solver. In particular, this numerical tool is applied to solve the fluid equation for a submerged body. The equations are solved for the following effects in the heave direction on a buoy: excitation force $F_e$, radiation force $F_r$ and restoring force $F_{res}$. Thus, the frequency domain modeling is performed with the Nemoh BEM solver. For each panel of a mesh (see Figure 2a), the hydrodynamic parameters are calculated for a frequency range. The pressures caused by the incoming wave are integrated into resulting forces on the buoy; see Figure 2b. From these data in the frequency domain, a filter is established to convert a wave into an exciting force (9). A magnitude diagram of the resulting filter is also shown in Figure 2b. The magnitude response resembles a second order filter with a cut-off frequency of $0.9$ rad/s.

An advanced function of Nemoh can be enabled to calculate the impulse response function (IRF). A comparison between the IRF identified and that provided by openWEC is shown in Figure 2c. To avoid convolution calculations in (2), it is replaced as ordinary differential equations (ODE) or as a transfer function in the Laplace domain (3). This can be performed using Prony’s method [21,31,32], which is implemented in MATLAB 2016b [33]. Different orders of ODEs have been tested to fit the impulse response. By comparing these responses, it was observed that for orders above eighth, the improvement in the system response was not significant.

$${\displaystyle \frac{F_e(s)}{\eta (s)}}={\displaystyle \frac{6.5\times {10}^5}{{(s+0.9)}^2}}\left[{\displaystyle \frac{N}{m}}\right]$$

(9)

$$K_r(s)={\displaystyle \frac{F_r(s)}{W(s)}}={\displaystyle \frac{a_{K_8}{s}^8+a_{K_7}{s}^7+a_{K_6}{s}^6+a_{K_5}{s}^5+a_{K_4}{s}^4+a_{K_3}{s}^3+a_{K_2}{s}^2+a_{K_1}s+a_{K_0}}{s^8+b_{K_7}{s}^7+b_{K_6}{s}^6+b_{K_5}{s}^5+b_{K_4}{s}^4+b_{K_3}{s}^3+b_{K_2}{s}^2+b_{K_1}s+b_{K_0}}}\left[{\displaystyle \frac{N}{m/s}}\right]$$

(10)

where the coefficients $a_K$ and $b_K$ are listed in

Table 1. Model coefficients identified for the WEC system.

Coefficient	Value	Coefficient	Value	Coefficient	Value	Coefficient	Value
—	—	$b_9$	8.8 × 101	—	—	—	—
$a_8$	2.4 × 10−6	$b_8$	1.5 × 105	$a_{K_8}$	1.2 × 103	—	—
$a_7$	2.2 × 10−4	$b_7$	6.7 × 10−6	$a_{K_7}$	2.3 × 105	$b_{K_7}$	8.9 × 101
$a_6$	3.6 × 10−1	$b_6$	4.5 × 109	$a_{K_6}$	1.9 × 108	$b_{K_6}$	1.5 × 105
$a_5$	1.6 × 101	$b_5$	1.3 × 1010	$a_{K_5}$	2.6 × 1010	$b_{K_5}$	6.7 × 106
$a_4$	1.1 × 104	$b_4$	6.4 × 1010	$a_{K_4}$	6.3 × 1012	$b_{K_4}$	4.5 × 109
$a_3$	3.3 × 104	$b_3$	8.6 × 1010	$a_{K_3}$	5.6 × 1014	$b_{K_3}$	1.3 × 1010
$a_2$	1.3 × 105	$b_2$	1.8 × 1011	$a_{K_2}$	1.7 × 1015	$b_{K_2}$	5.3 × 1010
$a_1$	1.4 × 105	$b_1$	1.1 × 1011	$a_{K_1}$	4.7 × 1015	$b_{K_1}$	5.5 × 1010
$a_0$	1.6 × 105	$b_0$	1.3 × 1011	$a_{K_0}$	1.4 × 1015	$b_{K_0}$	6.5 × 1010

.

Then, regrouping terms according to Equation (1), the transfer function (11) that relates external forces with the position of the buoy (z) is obtained.

$$G_{WEC}(s)={\displaystyle \frac{Z(s)}{F_{ext}(s)}}={\displaystyle \frac{a_8{s}^8+a_7{s}^7+a_6{s}^6+a_5{s}^5+a_4{s}^4+a_3{s}^3+a_2{s}^2+a_{1}s+a_0}{s^{10}+b_9{s}^9+b_8{s}^8+b_7{s}^7+b_6{s}^6+b_5{s}^5+b_4{s}^4+b_3{s}^3+b_2{s}^2+b_{1}s+b_0}}\left[{\displaystyle \frac{m}{N}}\right]$$

(11)

where the coefficients $a_n$ and $b_n$ are listed in Table 1.

2.3. Treatment of the Mathematical Model for the Design of MPCs

In order to ensure safe behavior and reduce mechanical fatigue in PAW systems, a model that allows them to constrain the oscillation speed and the position of the buoy is necessary. For this reason, a brief study of the identified model (11) is realized. By representing it in the state space, it can be verified that the system obtained is not completely controllable, and therefore, it is not a minimal realization [34]. Furthermore, it is not enough to reduce the model (11) to its minimum order. This is because, except for the system output (z), the other state variables lack physical sense, and this does not allow it to impose speed constraints (w) directly. As a solution, we study the transfer function (10) that defines the dynamics of the radiation force. Representing (10) in the state space, it can be seen how the model obtained is not of minimum order, so the system is minimized to get (12). Note that the state variables of this model ($x_r$) will not be controlled, so it is not necessary that they have physical sense.

$$\begin{array}{cc}\hfill {\dot{x}}_r(t)& =A_r{x}_r(t)+B_{r}w(t)\hfill \\ \hfill F{r}_{K_r}(t)& =C_r{x}_r(t)+D_{r}w(t)\hfill \end{array}$$

(12)

where the matrices $A_r$, $B_r$, $C_r$ and $D_r$ are given by (13).

$$\begin{array}{cc}\hfill A_r=& \left(\begin{array}{cccc}-3.0044& -1.4736& -0.3820& -0.2258\\ 8.1656& 0.0180& 0.0041& 0.0025\\ 0.0015& 4.0015& 0.0004& 0.0002\\ -0.0000& -0.0000& 2.0000& -0.0000\end{array}\right),B_r=\left(\begin{array}{c}121.3816\\ -1.3375\\ -0.1201\\ 0.0005\end{array}\right)\hfill \\ \hfill C_r=& {10}^3\left(\begin{array}{cccc}1.0083& 0.3683& 0.2624& 0.0377\end{array}\right),D_r=1218.70\hfill \end{array}$$

(13)

A complete model must consider the dynamics of the power take-off system. Given the similarity between the Wavestar system and the WEC studied in this work, the PTO model proposed in [24] is used, where the PTO dynamics are modeled according to (8). By representing this model in the state space, with the parameters indicated in [24], the matrices (14) have been obtained. For this purpose, the observable canonical form has been chosen. Thus, one of the state variables corresponds to the output of the PTO system, so MPCs may impose restrictions on the output of the actuator.

$$A_{pto}=\left(\begin{array}{cc}-8.7965& 1.0000\\ -157.9137& 0\end{array}\right),B_{pto}=\left(\begin{array}{c}0\\ 157.9137\end{array}\right),C_{pto}=\left(\begin{array}{cc}1& 0\end{array}\right),D_{pto}=0$$

(14)

After that, in this paper, we propose (15) as one of the design models, and the MPCs can impose restrictions on the following state variables: force applied by the PTO ($x_{pt{o}_1}$), oscillation speed (w) and buoy position (z). This model is similar to the one proposed in [8,12], but also considers the PTO dynamics.

$$\begin{array}{cc}\hfill \underset{\dot{x}}{\underbrace{\left(\begin{array}{c}\dot{z}\\ \dot{w}\\ {\dot{x}}_r\\ {\dot{x}}_{pto}\end{array}\right)}}& =\underset{A_{WEC}}{\underbrace{\left(\begin{array}{cccc}0& 1& 0& 0\\ -{\displaystyle \frac{k_{res}}{m_T}}& -{\displaystyle \frac{1}{m_T}}{D}_r& -{\displaystyle \frac{1}{m_T}}{C}_r& {\displaystyle \frac{1}{m_T}}{C}_{pto}\\ 0& B_r& A_r& 0\\ 0& 0& 0& A_{pto}\end{array}\right)}}\underset{x}{\underbrace{\left(\begin{array}{c}z\\ w\\ x_r\\ x_{pto}\end{array}\right)}}+\underset{B_{WE{C}_u}}{\underbrace{\left(\begin{array}{c}0\\ 0\\ 0\\ B_{pto}\end{array}\right)}}{F}_u+\underset{B_{WE{C}_{F_e}}}{\underbrace{\left(\begin{array}{c}0\\ {\displaystyle \frac{1}{m_T}}\\ 0\\ 0\end{array}\right)}}{F}_e\hfill \\ \hfill \underset{y}{\underbrace{\left(\begin{array}{c}z\\ w\end{array}\right)}}& =\underset{C_{WEC}}{\underbrace{\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right)}}\underset{x}{\underbrace{\left(\begin{array}{c}z\\ w\\ x_r\\ x_{pto}\end{array}\right)}}\hfill \end{array}$$

(15)

where $m_T=m+m_{\infty }$ and the matrices I (unit matrix) and zero have the required size according to their location, and the parameters not yet presented are listed in

Table 2. Parameters of the models (15) and (16) for the WEC system.

Symbol	Description	Value
$k_{res}$	Hydrostatic restoring coefficient	809,325 N/m
$B_{aprox}$	Stationary approximation to system damping	21,497 N s/m
m	Mass of water displaced by the buoy at rest	241,601.9 kg
$m_{\infty }$	Mass of water added	167,700 kg

.

Simplified Model for the Design

In this paper, as in others works [9,10,11,13,14], a simplified model (16) is used for the design of some MPCs. It differs from the previous model (15) in that it does not take into account the dynamic of the radiation force or the PTO dynamics. Figure 3 shows a comparison of the outputs of the simplified model and the complete model.

$$\begin{array}{cc}\hfill \underset{\dot{x}}{\underbrace{\left(\begin{array}{c}\dot{z}\\ \dot{w}\end{array}\right)}}& =\underset{A_{WECr}}{\underbrace{\left(\begin{array}{cc}0& 1\\ -{\displaystyle \frac{k_{res}}{m_T}}& -{\displaystyle \frac{B_{aprox}}{m_T}}\end{array}\right)}}\underset{x}{\underbrace{\left(\begin{array}{c}z\\ w\end{array}\right)}}+\underset{B_{WECr}}{\underbrace{\left(\begin{array}{c}0\\ {\displaystyle \frac{1}{m_T}}\end{array}\right)}}(F_u+F_e)\hfill \\ \hfill \underset{y}{\underbrace{\left(\begin{array}{c}z\\ w\end{array}\right)}}& =\underset{C_{WECr}}{\underbrace{\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)}}\underset{x}{\underbrace{\left(\begin{array}{c}z\\ w\end{array}\right)}}\hfill \end{array}$$

(16)

where $m_T=m+m_{\infty }$, and the parameter values are listed in Table 2.

3. Model Predictive Control for the Point Absorber WEC

This section details the design of the commonly-used MPCs for PAWs. Moreover, in this work, two MPCs are proposed based on the addition of an embedded integrator. In order to make a complete design, all the MPCs take into account constraints and the possibility of relaxing them, in the case of non-feasibility, in their cost functions. In particular, these constraints are applied to the control force of the PTO system, the position of the buoy and its oscillation speed. However, in the case of non-feasibility, only soft-constraints are applied to the position and oscillation speed of the system; due to the PTO being an actuator whose physical limit cannot be exceeded. Finally, for each controller, a sampling period of $T_m$ is set, which is used to discretize the mathematical model using the zero order hold (ZOH) approximation.

3.1. $MP{C}_1$

The cost function that minimizes this controller considers the maximization of extracted power directly; as the design model is used (16), for which the state vector estimated for a prediction horizon M and a control horizon N is defined according to Equation (17) [18,19].

$$X=\underset{J{x}_{(Mn\times n)}}{\underbrace{\left[\begin{array}{c}A\\ A^2\\ A^3\\ ⋮\\ A^M\end{array}\right]}}{x}_k+\underset{J{u}_{(Mn\times N)}}{\underbrace{\left[\begin{array}{ccccc}B& 0& 0& \dots & 0\\ AB& B& 0& \dots & 0\\ A^{2}B& AB& B& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A^{M-1}B& A^{M-2}B& A^{M-3}B& \dots & A^{M-N}B\end{array}\right]}}(F_{pto}+F_e)$$

(17)

where X represents the estimated state vector for a prediction horizon M, $x_k$ represents the state vector at the current instant, the matrices A and B are obtained from the model (16) discretized (ZOH), n is the order of the model and $F_{pto}$ and $F_e$ are vectors that contain the force applied by the PTO and the excitation force for whole control horizon N, respectively.

In a more compact form, the above equation can be expressed as:

$$X=J_x{x}_k+J_u(F_{pto}+F_e)$$

(18)

On the other side, as defined in [10,12,26], the mechanical power generated is given by (19). The expression of the power generated for the whole prediction horizon is obtained (20).

$$P_{gen}(t)=-w(t)F_{pto}(t)$$

(19)

$$P_{gen}=-W^T{F}_{pto}$$

(20)

where W and $F_{pto}$ are vectors with length M, which represent the oscillation speed and control force for the whole prediction horizon, respectively.

By replacing (18) in (20),

$$\begin{array}{cc}\hfill P_{gen}& =-{(S_{w}X)}^T{F}_{pto}\hfill \\ \hfill P_{gen}& =-{(S_w(J_x{x}_k+J_u{F}_{pto}+J_u{F}_e))}^T{F}_{pto}\hfill \end{array}$$

(21)

where $S_w$ is a selector matrix for speed w (size $M\times Mn$).

Developing (21) and grouping in terms of least squares, the cost function is obtained:

$$\begin{array}{cc}\hfill J(F_{pto})& ={\displaystyle \frac{1}{2}}{F_{pto}}^T\underset{H}{\underbrace{({J_u}^T{S_w}^T+R)}}{F}_{pto}+{\displaystyle \frac{1}{2}}\underset{b}{\underbrace{({F_e}^T{J_u}^T{S_w}^T+{x_k}^T{J_x}^T{S_w}^T)}}{F}_{pto}\hfill \end{array}$$

(22)

where the matrix R weights the control effort.

In addition, constraints are imposed to: force demanded for the PTO position and oscillation speed of the buoy (24). Therefore, the cost function (22) with constraints is defined as:

$$\begin{array}{cc}\hfill J(F_{pto})=& {\displaystyle \frac{1}{2}}{F_{pto}}^{T}H{F}_{pto}+{\displaystyle \frac{1}{2}}b{F}_{pto}\hfill \\ \hfill A_g{F}_{pto}\le & B_g\hfill \end{array}$$

(23)

where $A_g={\left[A_1{A}_2{A}_3\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3\right]}^T$; see Equation (24).

$$\begin{array}{cc}\hfill \underset{A_1}{\underbrace{\left[\begin{array}{c}I\\ -I\end{array}\right]}}{F}_{pto}\le \underset{B_1}{\underbrace{\left[\begin{array}{c}{F}_{PT{O}_{max}}\\ -F_{PT{O}_{min}}\end{array}\right]}},& \underset{A_2}{\underbrace{\left[\begin{array}{c}{S}_z{J}_u\\ -S_z{J}_u\end{array}\right]}}{F}_{pto}\le \underset{B_2}{\underbrace{\left[\begin{array}{c}{Z}_{n_{max}}-S_z(J_x{x}_k+J_u{F}_e)\\ -Z_{n_{min}}+S_z(J_x{x}_k+J_u{F}_e)\end{array}\right]}}\hfill \\ \hfill \underset{A_3}{\underbrace{\left[\begin{array}{c}{S}_w{J}_u\\ -S_w{J}_u\end{array}\right]}}{F}_{pto}\le & \underset{B_3}{\underbrace{\left[\begin{array}{c}{W}_{n_{max}}-S_w(J_x{x}_k+J_u{F}_e)\\ -W_{n_{min}}+S_w(J_x{x}_k+J_u{F}_e)\end{array}\right]}}\hfill \end{array}$$

(24)

where $S_z$ is a selector matrix for position (size $M\times Mn$), the vectors $F_{PT{O}_{max}}$ and $F_{PT{O}_{min}}$ (size $N\times 1$) define the nominal limits of the force applied by the PTO and the vectors $W_{n_{max}}$, $W_{n_{min}}$, $Z_{n_{max}}$ and $Z_{n_{min}}$ (size $M\times 1$) define the nominal limits of the buoy position and oscillation speed, respectively. In addition, the matrices I (unit matrix) and 0 have the required size according to their location.

In the case of non-feasibility, soft constraints to the position and oscillation speed of the system are applied. Thus, the cost function (22) would be defined as:

$$J(F_{pto},{\epsilon }_z,{\epsilon }_w)={\displaystyle \frac{1}{2}}{F_{pto}}^{T}H{F}_{pto}+b{F}_{pto}+{{\epsilon }_z}^T{W}_{{\epsilon }_z}{\epsilon }_z+{{\epsilon }_w}^T{W}_{{\epsilon }_w}{\epsilon }_w$$

(25)

where ${\epsilon }_z$ and ${\epsilon }_w$ represent the relaxation applied to position and speed along the prediction horizon M and the matrices $W_{{\epsilon }_z}$ and $W_{{\epsilon }_w}$ (size $M\times M$) weight these slacks.

By regrouping terms and adding soft restrictions, the cost function (25) can be expressed as:

$$\begin{array}{cc}\hfill J(\beta )=& {\displaystyle \frac{1}{2}}{\beta }^T\left[\begin{array}{ccc}H& 0& 0\\ 0& W_{{\epsilon }_z}& 0\\ 0& 0& W_{{\epsilon }_w}\end{array}\right]\beta +\left[\begin{array}{ccc}b& 0& 0\end{array}\right]\beta \hfill \\ \hfill A_g\beta \le & B_g\hfill \end{array}$$

(26)

where $\beta ={\left[F_{pto}{\epsilon }_z{\epsilon }_w\right]}^T$, with size $(N+2M)\times 1$, $A_g={\left[A_1{A}_2{A}_3{A}_4{A}_5\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3{B}_4{B}_5\right]}^T$; see Equation (27). Matrices 0 have the required size according to their location.

$$\begin{array}{c}\underset{A_2}{\underbrace{\left[\begin{array}{ccc}{S}_z{J}_u& -I& 0\\ -S_z{J}_u& -I& 0\end{array}\right]}}\beta \le \underset{B_2}{\underbrace{\left[\begin{array}{c}{Z}_{n_{max}}-S_z(J_x{x}_k+J_u{F}_e)\\ -Z_{n_{min}}+S_z(J_x{x}_k+J_u{F}_e)\end{array}\right]}},\underset{A_3}{\underbrace{\left[\begin{array}{ccc}{S}_w{J}_u& 0& I\\ -S_w{J}_u& 0& I\end{array}\right]}}\beta \le \underset{B_3}{\underbrace{\left[\begin{array}{c}{W}_{n_{max}}-S_w(J_x{x}_k+J_u{F}_e)\\ -W_{n_{min}}+S_w(J_x{x}_k+J_u{F}_e)\end{array}\right]}}\hfill \\ \underset{A_1}{\underbrace{\left[\begin{array}{ccc}I& 0& 0\\ -I& 0& 0\end{array}\right]}}\beta \le \underset{B_1}{\underbrace{\left[\begin{array}{c}{F}_{PT{O}_{max}}\\ -F_{PT{O}_{min}}\end{array}\right]}},\underset{A_4}{\underbrace{\left[\begin{array}{ccc}0& I& 0\\ 0& -I& 0\end{array}\right]}}\beta \le \underset{B_4}{\underbrace{\left[\begin{array}{c}{\kappa }_z\\ 0\end{array}\right]}},\underset{A_5}{\underbrace{\left[\begin{array}{ccc}0& 0& I\\ 0& 0& -I\end{array}\right]}}\beta \le \underset{B_5}{\underbrace{\left[\begin{array}{c}{\kappa }_w\\ 0\end{array}\right]}}\hfill \end{array}$$

(27)

where ${\kappa }_z$ and ${\kappa }_w$ are vectors (size $M\times 1$) that represent the maximum slack allowed for position and oscillation speed, respectively. The matrices I (unit matrix) and 0 have the required size according to their location.

3.2. $MP{C}_2$

This controller uses the simplified model (16). Its cost function is based on maximizing the extracted power by tracking a setpoint for the oscillation speed ($w_{ref}$) along a prediction horizon M,

$$J={(\tilde{w}-w_{ref})}^{T}Q(\tilde{w}-w_{ref})+{F_{pto}}^{T}R{F}_{pto}$$

(28)

where Q and R are diagonal matrices of size $(M\times M)$ and $(N\times N)$ that weight the tracking error and the control effort, respectively.

Substituting (18) in (28), the cost function for this controller can be written as (29).

$$J(F_{pto})={(\underset{f}{\underbrace{J_x{x}_k+J_u{F}_{pto}+J_u{F}_e}}-w_{ref})}^{T}Q(\underset{f}{\underbrace{J_x{x}_k+J_u{F}_{pto}+J_u{F}_e}}-w_{ref})+{F_{pto}}^{T}R{F}_{pto}$$

(29)

By developing the cost function (29) and grouping terms, the following expression is obtained:

$$J(F_{pto})={F_{pto}}^T\underset{H}{\underbrace{(J_u^T\delta J_u+R)}}{F}_{pto}+2\underset{b}{\underbrace{{(f-w_{ref})}^{T}Q{J}_u}}{F}_{pto}+\underset{l}{\underbrace{{(f-w_{ref})}^{T}Q(f-w_{ref})}}$$

(30)

Note that the term l can be ignored when the cost function is minimized, because it does not depend on the variable to be optimized ($F_{pto}$),

$$J(F_{pto})={\displaystyle \frac{1}{2}}{F_{pto}}^T\underset{H}{\underbrace{(J_u^T\delta J_u+R)}}{F}_{pto}+\underset{b}{\underbrace{{(J_x{x}_k+J_u{F}_{pto}+J_u{F}_e-w_{ref})}^{T}Q{J}_u}}{F}_{pto}$$

(31)

The constraints imposed on this cost function can be expressed in the same way as in the $MP{C}_1$ controller; using (23) for hard constraints and (26) for soft constraints. On the other hand, the reference trajectory, or setpoint for the oscillation speed, is defined by the approach proposed in [4], $w_{ref}=F_e/2B_{aprox}$.

3.3. $MP{C}_3$

This approach is a contribution made in this work. This controller uses the simplified model (16) to which an embedded integrator has been added according to the theory of predictive controllers in the state space [18,19]. Its cost function is based on the maximization of the extracted power through the tracking of a setpoint for the oscillation speed. To add the integrator, it is necessary to multiply the model (16) discretized by the operator $▵=1-z^{-1}$. Regrouping terms, an extended state vector is defined as:

$$\begin{array}{cc}\hfill \underset{x_e(t+1)}{\underbrace{\left[\begin{array}{c}▵x(t+1)\\ y(t+1)\end{array}\right]}}& =\underset{A_e}{\underbrace{\left[\begin{array}{cc}A& 0\\ CA& I\end{array}\right]}}\underset{x_e(t)}{\underbrace{\left[\begin{array}{c}▵x(t)\\ y(t)\end{array}\right]}}+\underset{B_e}{\underbrace{\left[\begin{array}{c}B\\ CB\end{array}\right]}}\left(▵F_{pto}(t)+▵F_e(t)\right)\hfill \\ \hfill y(t)& =\underset{C_e}{\underbrace{\left[\begin{array}{cc}0& I\end{array}\right]}}\underset{x_e(t)}{\underbrace{\left[\begin{array}{c}▵x(t)\\ y(t)\end{array}\right]}}\hfill \end{array}$$

(32)

where the output vector $y(t)$ is formed by the position and oscillation speed of the WEC system, $▵x(t)$ represents the state vector increment, the matrices A, B and C are from the model (16) discretized (ZOH) and the matrices 0 have the required size according to their location.

Using the extended model (32), the prediction of the outputs (z, w) is defined for a prediction horizon M and a control horizon N according to:

$$Y=\underset{F_{(2M\times n_e})}{\underbrace{\left[\begin{array}{c}CA\\ C{A}^2\\ C{A}^3\\ ⋮\\ C{A}^M\end{array}\right]}}{X}_e+\underset{G_{(2M\times N})}{\underbrace{\left[\begin{array}{ccccc}CB& 0& 0& \dots & 0\\ CAB& CB& 0& \dots & 0\\ C{A}^{2}B& CAB& CB& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C{A}^{M-1}B& C{A}^{M-2}B& C{A}^{M-3}B& \dots & C{A}^{M-N}B\end{array}\right]}}(▵F_{pto}+▵F_e)$$

(33)

where $n_e=n+j$ represents the order of the extended model and j its number of outputs. In the matrices A, B and C, the sub-index e has been omitted to get a clearer notation.

By adding the embedded integrator, this controller minimizes a cost function that gets the optimal increase in control force ($F_{pto}$) for the full control horizon N,

$$J={(\tilde{w}-w_{ref})}^{T}Q(\tilde{w}-w_{ref})+{▵F_{pto}}^{T}R▵F_{pto}$$

(34)

where Q and R are diagonal matrices of size $(M\times M)$ and $(N\times N)$ that weigh the tracking error and the control effort, respectively.

Replacing the output prediction (33) in (34) and obviating the independent term of $▵F_{pto}$:

$$J(▵F_{pto})={\displaystyle \frac{1}{2}}{▵F_{pto}}^T\underset{H}{\underbrace{(G^T\delta G+R)}}▵F_{pto}+\underset{b}{\underbrace{{(F{x}_{e_k}+G▵F_e-w_{ref})}^{T}QG}}▵F_{pto}$$

(35)

In addition, constraints are added to the demanded force on the PTO, position and oscillation speed of the system. Therefore, the cost function (35) subject to the restrictions is defined as:

$$\begin{array}{cc}\hfill J(▵F_{pto})=& {\displaystyle \frac{1}{2}}{▵F_{pto}}^{T}H▵F_{pto}+b▵F_{pto}\hfill \\ \hfill A_g▵F_{pto}\le & B_g\hfill \end{array}$$

(36)

where $A_g={\left[A_1{A}_2{A}_3\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3\right]}^T$; see Equation (37).

$$\begin{array}{c}\underset{A_1}{\underbrace{\left[\begin{array}{c}T\\ -T\end{array}\right]}}▵F_{pto}\le \underset{B_1}{\underbrace{\left[\begin{array}{c}{F}_{PT{O}_{max}}\\ -F_{PT{O}_{min}}\end{array}\right]}},\underset{A_2}{\underbrace{\left[\begin{array}{c}{S}_{z}G\\ -S_{z}G\end{array}\right]}}▵F_{pto}\le \underset{B_2}{\underbrace{\left[\begin{array}{c}{Z}_{n_{max}}-S_z(F{x}_{e_k}+G▵F_e)\\ -Z_{n_{min}}+S_z(F{x}_{e_k}+G▵F_e)\end{array}\right]}}\\ \underset{A_3}{\underbrace{\left[\begin{array}{c}{S}_{w}G\\ -S_{w}G\end{array}\right]}}▵F_{pto}\le \underset{B_3}{\underbrace{\left[\begin{array}{c}{W}_{n_{max}}-S_w(F{x}_{e_k}+G▵F_e)\\ -W_{n_{min}}+S_w(F{x}_{e_k}+G▵F_e)\end{array}\right]}}\end{array}$$

(37)

where $S_z$ and $S_w$ are selector matrices for z and w (size $M\times Mn$), T is a lower triangular matrix (size $M\times N$), the vectors $F_{PT{O}_{max}}$ and $F_{PT{O}_{min}}$ (size $N\times 1$) define the nominal limits of the force applied by the PTO and the vectors $W_{n_{max}}$, $W_{n_{min}}$, $Z_{n_{max}}$ and $Z_{n_{min}}$ (size $M\times 1$) define the nominal limits of the buoy position and oscillation speed, respectively. The matrices I (unit matrix) and 0 have the required size according to their location.

In the case of non-feasibility in the cost function (36), soft constraints are applied,

$$\begin{array}{cc}\hfill J(\beta )=& {\displaystyle \frac{1}{2}}{\beta }^T\left[\begin{array}{ccc}H& 0& 0\\ 0& W_{{\epsilon }_z}& 0\\ 0& 0& W_{{\epsilon }_w}\end{array}\right]\beta +\left[\begin{array}{ccc}b& 0& 0\end{array}\right]\beta \hfill \\ \hfill A_g\beta \le & B_g\hfill \end{array}$$

(38)

where $A_g={\left[A_1{A}_2{A}_3{A}_4{A}_5\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3{B}_4{B}_5\right]}^T$; see Equation (38).

$$\begin{array}{c}\underset{A_1}{\underbrace{\left[\begin{array}{ccc}T& 0& 0\\ -T& 0& 0\end{array}\right]}}\beta \le \underset{B_1}{\underbrace{\left[\begin{array}{c}{F}_{PT{O}_{max}}\\ -F_{PT{O}_{min}}\end{array}\right]}}\hfill \\ \underset{A_2}{\underbrace{\left[\begin{array}{ccc}{S}_{z}G& -I& 0\\ -S_{z}G& -I& 0\end{array}\right]}}\beta \le \underset{B_2}{\underbrace{\left[\begin{array}{c}{Z}_{n_{max}}-S_z(F{x}_{e_k}+G▵F_e)\\ -Z_{n_{min}}+S_z(F{x}_{e_k}+G▵F_e)\end{array}\right]}},\underset{A_4}{\underbrace{\left[\begin{array}{ccc}0& I& 0\\ 0& -I& 0\end{array}\right]}}\beta \le \underset{B_4}{\underbrace{\left[\begin{array}{c}{\kappa }_z\\ 0\end{array}\right]}}\hfill \\ \underset{A_3}{\underbrace{\left[\begin{array}{ccc}{S}_{w}G& 0& I\\ -S_{w}G& 0& I\end{array}\right]}}\beta \le \underset{B_3}{\underbrace{\left[\begin{array}{c}{W}_{n_{max}}-S_w(F{x}_{e_k}+G▵F_e)\\ -W_{n_{min}}+S_w(F{x}_{e_k}+G▵F_e)\end{array}\right]}},\underset{A_5}{\underbrace{\left[\begin{array}{ccc}0& 0& I\\ 0& 0& -I\end{array}\right]}}\beta \le \underset{B_5}{\underbrace{\left[\begin{array}{c}{\kappa }_w\\ 0\end{array}\right]}}\hfill \end{array}$$

(39)

where ${\kappa }_z$ and ${\kappa }_w$ are vectors (size $M\times 1$) that represent the maximum slack allowed for z and w, respectively. The matrices I (unit matrix) and 0 have the required size according to their location.

3.4. $MP{C}_4$

This controller is made using the model (15). The cost function that minimizes this controller is focused on the maximization of extracted power directly. The matrix development needed to express this controller as a least squares problem is analogous to that presented for controller $MP{C}_1$. Although, when the PTO dynamics are taken into account, Equation (18) should be redefined as:

$$X=J_x{x}_k+J_u{F}_{pto}+J_f{F}_e$$

(40)

where $J_x$ is a matrix already defined in Equation (17), while the matrices $J_u$ and $J_f$ are given by:

$$\begin{array}{cc}\hfill J_u& =\left[\begin{array}{ccccc}{B}_u& 0& 0& \dots & 0\\ A{B}_u& B_u& 0& \dots & 0\\ A^2{B}_u& A{B}_u& B_u& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A^{M-1}{B}_u& A^{M-2}{B}_u& A^{M-3}{B}_u& \dots & A^{M-N}{B}_u\end{array}\right]\hfill \\ \hfill J_f& =\left[\begin{array}{ccccc}{B}_{F_e}& 0& 0& \dots & 0\\ A{B}_{F_e}& B_{F_e}& 0& \dots & 0\\ A^2{B}_{F_e}& A{B}_{F_e}& B_{F_e}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A^{M-1}{B}_{F_e}& A^{M-2}{B}_{F_e}& A^{M-3}{B}_{F_e}& \dots & A^{M-N}{B}_{F_e}\end{array}\right]\hfill \end{array}$$

(41)

where A, $B_u$ and $B_{F_e}$ are obtained by discretizing (ZOH) $A_{WEC}$, $B_{WE{C}_u}$ and $B_{WE{C}_{F_e}}$ of the model (15).

The cost function to be minimized by this controller can be expressed according to (42). Its development is analogous to that carried out for the controller $MP{C}_1$.

$$\begin{array}{cc}\hfill J(F_{pto})& ={\displaystyle \frac{1}{2}}{F_{pto}}^T\underset{H}{\underbrace{({J_u}^T{S_w}^T+R)}}{F}_{pto}+{\displaystyle \frac{1}{2}}\underset{b}{\underbrace{({F_e}^T{J_f}^T{S_w}^T+{x_k}^T{J_x}^T{S_w}^T)}}{F}_{pto}\hfill \end{array}$$

(42)

In addition, the nominal constraints imposed on the system must be added, which are defined in the same way as in the controller $MP{C}_1$, Equation (24). However, in this case, it must be considered that the state vector prediction is given by Equation (40). Therefore, the cost function (42) subject to the constraints is defined as:

$$\begin{array}{cc}\hfill J(F_{pto})=& {\displaystyle \frac{1}{2}}{F_{pto}}^{T}H{F}_{pto}+{\displaystyle \frac{1}{2}}b{F}_{pto}\hfill \\ \hfill A_g{F}_{pto}\le & B_g\hfill \end{array}$$

(43)

where $A_g={\left[A_1{A}_2{A}_3\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3\right]}^T$; see Equation (24).

Finally, in the case of non-feasibility in the function (48), soft constraints will be applied to the system. These can be expressed in a similar way to the development shown for $MP{C}_1$, Equation (27). However, in this case, it must be considered that the prediction of the state vector is given by Equation (40). Therefore, the cost function of this controller subject to soft constraints is given by:

$$\begin{array}{cc}\hfill J(\beta )=& {\displaystyle \frac{1}{2}}{\beta }^T\left[\begin{array}{ccc}H& 0& 0\\ 0& W_{{\epsilon }_z}& 0\\ 0& 0& W_{{\epsilon }_w}\end{array}\right]\beta +\left[\begin{array}{ccc}b& 0& 0\end{array}\right]\beta \hfill \\ \hfill A_g\beta \le & B_g\hfill \end{array}$$

(44)

where $A_g={\left[A_1{A}_2{A}_3{A}_4{A}_5\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3{B}_4{B}_5\right]}^T$; see Equation (27).

3.5. MPC₅

This controller is another contribution of this work. It uses the model (15), to which an embedded integrator has been added according to the theory of predictive controllers in the state space [18,19]. Its cost function to minimize is based on the maximization of the extracted power through the tracking of a setpoint for the oscillation speed of the system $w_{ref}$. As in the $MP{C}_3$, the state vector is extended by adding an embedded integrator (32). Despite taking into account the PTO dynamics, the prediction of the output vector for the full prediction horizon M is defined by:

$$Y=G_u▵F_{pto}+\underset{f}{\underbrace{F{x}_{e_k}+G_f▵F_e}}$$

(45)

where F is a matrix already defined in Equation (33), and the matrices $G_u$ and $G_f$ are given by:

$$\begin{array}{cc}\hfill G_u& =\left[\begin{array}{ccccc}{C}_e{B}_{e_u}& 0& 0& \dots & 0\\ C_e{A}_e{B}_{e_u}& C_e{B}_{e_u}& 0& \dots & 0\\ C_e{A_e}^2{B}_{e_u}& C_e{A}_e{B}_{e_u}& C_e{B}_{e_u}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C_e{A_e}^{M-1}{B}_{e_u}& C_e{A_e}^{M-2}{B}_{e_u}& C_e{A_e}^{M-3}{B}_{e_u}& \dots & C_e{A_e}^{M-N}{B}_{e_u}\end{array}\right]\hfill \\ \hfill G_f& =\left[\begin{array}{ccccc}{C}_e{B}_{e_{F_e}}& 0& 0& \dots & 0\\ C_e{A}_e{B}_{e_{F_e}}& C_e{B}_{e_{F_e}}& 0& \dots & 0\\ C_e{A_e}^2{B}_{e_{F_e}}& C_e{A}_e{B}_{e_{F_e}}& C_e{B}_{e_{F_e}}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C_e{A_e}^{M-1}{B}_{e_{F_e}}& C_e{A_e}^{M-2}{B}_{e_{F_e}}& C_e{A_e}^{M-3}{B}_{e_{F_e}}& \dots & C_e{A_e}^{M-N}{B}_{e_{F_e}}\end{array}\right]\hfill \end{array}$$

(46)

where $B_{e_u}$ and $B_{e_{F_e}}$ are obtained by discretizing the matrices that define the inputs of (15) extended.

Analogous to controller $MP{C}_3$, the cost function to be minimized can be expressed according to:

$$J(▵F_{pto})={\displaystyle \frac{1}{2}}{▵F_{pto}}^T\underset{H}{\underbrace{({G_u}^T\delta G_u+R)}}▵F_{pto}+\underset{b}{\underbrace{{(F{x}_{e_k}+G_f▵F_e-w_{ref})}^{T}Q{G}_u}}▵F_{pto}$$

(47)

In addition, it is necessary to add nominal constraints to the system, which are defined as in the controller $MP{C}_3$. However, in this case, it must be considered that the prediction of the system outputs is given by (45). Thus, the cost function (47) subject to the constraint is defined as:

$$\begin{array}{cc}\hfill J(▵F_{pto})=& {\displaystyle \frac{1}{2}}{▵F_{pto}}^{T}H▵F_{pto}+{\displaystyle \frac{1}{2}}b▵F_{pto}\hfill \\ \hfill A_g▵F_{pto}\le & B_g\hfill \end{array}$$

(48)

where $A_g={\left[A_1{A}_2{A}_3\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3\right]}^T$; see Equation (37).

Furthermore, in the case of the non-feasibility in the function (47), soft constraints are used. The soft constraints are expressed analogously to the development made for controller $MP{C}_3$, Equation (39). However, in this case, the prediction of the system output is given by (45). Therefore, the cost function of this controller subject to soft constraints is defined as:

$$\begin{array}{cc}\hfill J(\beta )=& {\displaystyle \frac{1}{2}}{\beta }^T\left[\begin{array}{ccc}H& 0& 0\\ 0& W_{{\epsilon }_z}& 0\\ 0& 0& W_{{\epsilon }_w}\end{array}\right]\beta +\left[\begin{array}{ccc}b& 0& 0\end{array}\right]\beta \hfill \\ \hfill A_g\beta \le & B_g\hfill \end{array}$$

(49)

where $A_g={\left[A_1{A}_2{A}_3{A}_4{A}_5\right]}^T$ and $B_g={\left[B_1{B}_2{B}_3{B}_4{B}_5\right]}^T$; see Equation (39).

3.6. Conventional Controllers

In order to make a comparative analysis of the proposed predictive controllers, this section presents the design of two of the controllers most commonly used in WEC systems [5,6]. Firstly, a resistive damping (RD) controller (or proportional control) has been tuned for the PTO system, whose control law is given by:

$$F_{pto}(k)=-K_{RD}w(k)$$

(50)

where k is the current sampling time and $K_{RD}$ is a proportional gain ($K_P$), which must be tuned to maximize the power generated.

On the other hand, an I-P control has been designed whose control law is defined by Equation (51). A compromise between generated power and keeping the system within its nominal limits of operation is maintained by tuning the gains $K_P$ and $K_I$.

$$\begin{array}{cc}\hfill F_{pto}(k)& =u_P(k)+u_I(k)\hfill \\ \hfill u_P(k)& =-K_{P}w(k)\hfill \\ \hfill u_I(k)& =K_I{T}_m(w_{ref}(k)-w(k))+u_I(k-1)\hfill \end{array}$$

(51)

where k is the current moment. The backward Euler method has been used to discretize the controller for a sampling period $T_m$.

4. Study about Performances and Robustness

This section presents the results obtained from an in-depth study of the performance and robustness of the five MPCs designs. This study assumes that the system state vector and the excitation force ($F_e$) for the whole prediction horizon (M) are known. The computer simulations have been carried out using Simulink, employing a fourth order Runge–Kutta integration method with a fixed step of one millisecond. In order to solve the optimization problems with constraints, associated with MPCs’ design, the MATLAB quadprog function has been employed [35]. Note that although the quadprog function provides $F_{pto}$ for the whole control horizon N, only the setpoint obtained for the current instant k is applied [18,19].

Due to the fact that the WEC system of the W2POWER platform has not yet been built, its operating limits and physical limits are not available. Therefore, these limits have been chosen on the basis of [12,15,16,24], whose WEC systems are similar to the system studied in this paper; see

Table 3. Hard and soft constraints for the WEC system.

Symbol	Description	Value
$F_{PT{O}_{max}}$	Maximum stationary force for the power take-off system	450 $KN$
$F_{PT{O}_{min}}$	Minimum stationary force for the power take-off system	−450 $KN$
$z_{n_{max}}$	Maximum nominal limit for the buoy position	1.25 m
$z_{n_{min}}$	Minimum nominal limit for the buoy position	−1.25 m
$w_{n_{max}}$	Maximum nominal limit for oscillation speed	1 m/s
$w_{n_{min}}$	Minimum nominal limit for oscillation speed	−1 m/s
$z_{f_{max}}$	Maximum physical limit for the buoy position	1.7 m
$z_{f_{min}}$	Minimum physical limit for the buoy position	−1.7 m
$w_{f_{max}}$	Maximum physical limit for oscillation speed	1.3 m/s
$w_{f_{min}}$	Minimum physical limit for oscillation speed	−1.3 m/s
${\kappa }_z$	Maximum slack applied to the position nominal limit	0.45 m
${\kappa }_w$	Maximum slack applied to the speed nominal limit	0.3 m/s

.

4.1. Performance Comparison

This section shows a comparison between the seven controllers designed. The controllers’ performances are compared in terms of: average power generated, reduction of instantaneous power peaks, overshoot of nominal limits and control effort. A sea state defined by the JONSWAP spectrum [15,16], with a significant wave height of three meters and a peak period of 11 s, has been chosen as a realistic scenario for evaluating these characteristics. In order to obtain a truthful performance comparison, all predictive controllers must have the same information about the incoming wave, which is recorded in the prediction time $t_f$. To set the prediction horizon, two factors have been taken into account. Firstly, in [7,8,10], the authors used realistic sea states, and they set the prediction times between two and four seconds. Furthermore, in [8], the authors checked that $t_f$ could also be reduced to three or four seconds without significant reduction of the harvested energy. Secondly, in [36], the authors showed how short-term wave forecasting models maintain good performance up to five seconds of prediction for wide wave spectra. Therefore, in this work, the prediction time chosen was three seconds, the same as the one set in [10]. With this $t_f$, a prediction horizon $M=75$, a control horizon $N=75$ and a sampling period $T_m=0.04$ s has been set for all MPCs. The RD and I-P controllers have the same sampling period. On the other hand, as a result of a fine-tuning for this realistic sea state, the parameters of the controllers are listed in

Table 4. Control parameter set for the sea state defined by the JONSWAP spectrum (3 m of significant wave height and 11 s of peak period).

Controller	R	Q	${\mathit{W}}_{{\mathit{\epsilon }}_{\mathit{z}}}$	${\mathit{W}}_{{\mathit{\epsilon }}_{\mathit{w}}}$	${\mathit{K}}_{\mathit{P}}$	${\mathit{K}}_{\mathit{I}}$
MPC1	4.500 × 10−7	—	1.000 × 1010	1.000 × 107	—	—
MPC2	1.000 × 10−7	9.150 × 104	1.000 × 107	2.000 × 109	—	—
MPC3	5.000 × 105	4.575 × 104	5.000 × 1010	1.000 × 105	—	—
MPC4	4.500 × 10−7	—	1.000 × 108	1.000 × 104	—	—
MPC5	5.000 × 10−5	4.000 × 105	1.000 × 107	1.000 × 109	—	—
RD	—	—	—	—	7.902 × 105	—
I-P	—	—	—	—	7.050 × 105	2.228 × 104

.

Figure 4 shows a comparison of the instantaneous mechanical powers generated by applying the seven controllers to the mathematical model (15). In this comparison, it can be seen how MPCs with an embedded integrator (MPC₃ and MPC₅) achieve more regular power than the MPCs most used for WEC systems, those whose optimization criteria maximize the extracted power directly (MPC₁ and MPC₄). In addition, Figure 5 shows how MPCs with embedded integrators achieve less overshoot in the control force applied by the PTO system than all other MPCs (reducing actuator overstress).

By applying the MPC₂ to the system, it gives the most irregular power; while the power generated with the MPC₃, designed from the same model and following the same optimization criteria, is much cleaner, with fewer occasional peaks. This is due to the fact that the MPC₁ is continuously applying soft constraints to the oscillation speed, because it does not carry out a good control of the force that the PTO system exerts on the WEC during the time (see Figure 5). On the other hand,

Table 5. Results obtained in the application of the seven controllers to the WEC system. The powers are expressed in kW. Note that $ONLP$ indicates overshoot of nominal limits for position and $ONLS$ indicates overshoot of nominal limits for speed.

Controller	${\overline{\mathit{P}}}_{\mathit{gen}}$	$\mathit{ONLP}$	$\mathit{ONLS}$	${\mathit{P}}_{{\mathit{gen}}_{\mathit{Max}}}$	${\mathit{P}}_{{\mathit{gen}}_{\mathit{Min}}}$
MPC1	127.60	0.0000	0.0023	437.48	−356.33
MPC2	110.72	0.0000	0.0281	744.06	−364.97
MPC3	125.53	0.0000	0.0000	444.24	−183.62
MPC4	127.50	0.0000	0.0107	452.75	−369.27
MPC5	129.01	0.0000	0.0413	481.53	−182.85
RD	67.02	0.0000	0.0000	258.88	−1.05
I-P	109.90	0.0073	0.0511	583.73	−107.34

shows how the addition of the embedded integrator (MPC₃) considerably improves the behavior of the system with respect to that obtained by applying the MPC₂. Since, the MPC₃ does not exceed nominal limits of the position and oscillation speed on any occasion. Furthermore, the MPC₃ controller generates a clearer control signal than the MPC₂ (see Figure 5), and as a consequence, the underdamped response of the PTO system decreases greatly.

Table 5 records the most significant quality indicators of the control performed by each controller. First, this table shows the average mechanical powers generated by the WEC system when applying each controller. In this aspect, the MPC₅ controller is the one that generates more power, followed very closely by the MPC₁, MPC₄ and, with a bit more distance, the MPC₃. On the other side, the indicators $ONLP$ and $ONLS$ quantify the area of overshoot from the nominal limits of the position and oscillation speed, respectively. In this sense, the MPC₃ (together with the RD control) provides the best behavior, since it does not apply slack to the nominal limits on any occasion, then it is followed by MPC₁. This can be verified in Figure 6 and Figure 7, which show a comparison between the positions and oscillation speeds obtained by applying the designed controllers to the mathematical model (15). As can be seen, all controllers keep the WEC system within its physical operating limits. Finally, the last two columns of Table 5 show the maximum and minimum power peaks obtained when applying each controller. In this aspect, ignoring the resistive damping control, MPCs with embedded integrators are once again the best performers. Note that these power peaks will cause an oversizing of: electrical machines, power electronics, accumulators, etc.

With respect to the two optimization criteria compared in this paper, it can be concluded that the MPCs that employ an optimization criterion based on the minimization of the error between w and $w_{ref}$ should only be used if they are designed from a model with an embedded integrator. Thus, this approach achieves better performance (in terms of diminution of instantaneous power peaks and reduction of mechanical fatigue due to exceeding nominal limits) than the standard optimization criteria based on maximizing the extracted power directly. On the other hand, with respect to the use of a complete model or a simplified model, it can be concluded that (in terms of instantaneous power generated and reduction of mechanical fatigue) the use of a complete model does not improve the performance of MPCs for this WEC significantly. This is because, the increase in the average power generated is minimal, and only in the case of the MPC₁ and MPC₄, the use of a complete model reduces the overshoot of the PTO system a bit. Finally, the variation of the extracted mechanical power as a function of the prediction time $t_f$ is studied. In particular, Figure 8 shows a comparison between MPC₃ and MPC₁. As can be seen, in both cases, after a prediction time of 3 s, the power generation does not improve significantly, while the computation effort increases. In addition, it should be pointed that MPC₃ does not have a monotonously increasing behavior that relates $t_f$ to the power generated, as would be expected. Nevertheless, in this aspect, the MPC₃ is better than the MPC₁, because it can generate more power with less information of the future excitation force (up to $2.5$ s).

To conclude the performance study, by comparing the previous figures and the values recorded in Table 5, it can be seen how the MPCs almost double the average power generated by the WEC system with respect to that obtained by tuning an adequate resistive damping for the PTO system (RD control). In addition, like the RD control, the MPC₃ keeps the system within its nominal operating limits. On the other hand, it is also verified that the MPCs are superior in performances than the conventional I-P controller, in terms of: average mechanical powers generated, diminution of instantaneous power peaks and reduction of mechanical fatigue (due to exceeding nominal limits).

4.2. Robustness Comparison

Another contribution of this work, searching for the greatest realism in the comparative of the designed controllers, is that uncertainty is added to the complete system model (15) to the most significant identified parameters (52): added mass and dynamics of the radiation force, which have been obtained through the openWEC software, and the hydrostatic restoring coefficient $K_{res}$ (a nonlinear parameter that has been linearized during the modeling of the WEC system).

$$F_{res}(t)=-k_{res}(1+{▵}_{k_{res}})z(t),F{r}_{m_{\infty }}(t)=m_{\infty }(1+{▵}_{m_{\infty }})\dot{w}(t),{\displaystyle \frac{F{r}_{K_r}(s)}{W(s)}}=(1+{▵}_B)K_r(s)$$

(52)

where Δ represents the added uncertainty in each parameter.

Note that when modifying the physical parameters of the WEC system, the frequency response of the filter (9) is affected. Therefore, looking for a truthful comparison, it would be necessary to identify a new filter (6) for each added uncertainty. Given the high number of simulations required, applying different levels of uncertainty to each of the parameters, this is not feasible. For this reason, in this work, the Morison model (7) is used to define $F_e$, allowing one to modify the excitation force caused by the wave as a function of the added uncertainty more easily. Thus, when modifying the physical parameters of the system, the external force that the wave causes on the system also varies. Therefore, Equation (7) is redefined for this analysis as:

$$F_e(t)=m_{\infty }(1+{▵}_{m_{\infty }})\ddot{\eta }(t)+B(1+{▵}_B)\dot{\eta }(t)+k_{res}(1+{▵}_{k_{res}})\eta (t)$$

(53)

It should be noted that Equation (53) uses a first and a second derivative of wave height. As a consequence, if the Morison model is used for the excitation force in a realistic sea state, it will amplify the high-frequency harmonics that are part of the wave. This is the opposite of the bandwidth of the WEC system provided by the openWEC software (see Figure 2b). For this reason, in this part of the paper, the simulations are performed in an irregular sea state formed by the fifteen sinusoidal components listed in

Table 6. Sinusoidal components used for sea-state (values expressed in international units).

Component	Amplitude	Period	Phase	Component	Amplitude	Period	Phase
$s_1$	$0.420$	$13.00$	$-\pi$	$s_9$	$0.200$	$9.00$	$0.00$
$s_2$	$0.520$	$12.50$	$1.5\pi$	$s_{10}$	$0.180$	$8.50$	$\pi$
$s_3$	$0.420$	$12.25$	$0.40$	$s_{11}$	$0.200$	$7.50$	$0.10$
$s_4$	$0.520$	$11.50$	$0.20$	$s_{12}$	$0.150$	$6.50$	$-0.77$
$s_5$	$0.450$	$11.25$	$0.11\pi$	$s_{13}$	$0.100$	$5.50$	$0.5\pi$
$s_6$	$0.300$	$10.50$	$-1.50$	$s_{14}$	$0.075$	$5.00$	$0.00$
$s_7$	$0.500$	$10.00$	$-0.33$	$s_{15}$	$0.020$	$3.70$	$0.12$
$s_8$	$0.210$	$9.50$	$0.78$	—	—	—	—

.

In order to create an unfavorable test scenario for the MPCs, two factors have been adjusted. In the first place, the sea state recorded in Table 6 is not favorable to the controllers, because the wave force becomes more than twice the stationary force that the PTO system can apply (recorded in Table 3). Moreover, the prediction time $t_f$ has been limited to a maximum of $1.5$ s. After this, a fine-tuning has been made to all the MPCs looking for a balance between the generated power and the overshoot of the nominal limits of the WEC system. The result of this fine-tuning is listed in

Table 7. Control parameters set for the sea state recorded in Table 7.

Controller	${\mathit{T}}_{\mathit{m}}$$\left[\mathit{s}\right]$	${\mathit{T}}_{\mathit{f}}\left[\mathit{s}\right]$	R	Q	${\mathit{W}}_{{\mathit{\epsilon }}_{\mathit{z}}}$	${\mathit{W}}_{{\mathit{\epsilon }}_{\mathit{w}}}$
MPC1	$0.05$	$1.5$	$1.1\times {10}^{-7}$	—	$1.0\times {10}^{10}$	$1.0\times {10}^7$
MPC2	$0.04$	$1.2$	$1.0\times {10}^{-7}$	$2.15\times {10}^4$	$1.0\times {10}^7$	$2.0\times {10}^9$
MPC3	$0.04$	$1.2$	$5.0\times {10}^{-5}$	$1.475\times {10}^4$	$5.0\times {10}^{10}$	$1.0\times {10}^5$
MPC4	$0.05$	$1.5$	$4.5\times {10}^{-7}$	—	$1.0\times {10}^8$	$1.0\times {10}^4$
MPC5	$0.02$	$0.6$	$5.0\times {10}^{-5}$	$1.75\times {10}^4$	$1.0\times {10}^7$	$1.0\times {10}^9$

.

Note that for the same wave height, the excitation force can increase or decrease according to the uncertainty added in each parameter. Therefore, there will be situations where, for the sea state defined in Table 6, the controller cannot keep the system within its physical limits. This is because, the actuation force of the PTO system will be much lower than the excitation force. In this paper, this non-feasibility situation will be considered as the robustness limit that the controller can support. This limit is defined for the uncertainty added in each of the parameters (52). This non-feasibility situation with soft constraints does not mean that the closed-loop system becomes unstable, but that the controller cannot keep the WEC system within the physical limits defined in Table 3. Furthermore, in order to obtain a more complete analysis of how this uncertainty affects the closed-loop system, the average powers generated for each value of added uncertainty to each parameter are recorded. A large number of simulations has been carried out for this purpose; all of them have a duration of 120 s and use fourth order Runge–Kutta (RK4) integration method with an integration step of 1 ms.

Once the robustness study has been defined, Figure 9 shows the results obtained as a function of the added uncertainties; feasible limits obtained for the five MPCs and the variation of their mean power generated. With respect to the added uncertainty in $m_{\infty }$, the MPC₂ is the least robust. Meanwhile, the MPC₁ offers the best features in a power-robustness ratio. However, it should be noted that when considering more reasonable added uncertainty values (interval $[-50,50]\%$), the MPC₅ extracts significantly more power than the others. It should also be noted that the controllers that directly maximize power in their cost function (MPC₁ and MPC₄) have the most predictable behavior with respect to the added uncertainty in $m_{\infty }$. On the other hand, the MPC₅ offers the best features with respect to the uncertainty added to the dynamics of the radiation force (up to $400\%$). In contrast, the MPC₂ gives very bad results in this respect. The MPC₁ also gets good results, because it achieves a practically constant power production despite variations of ${▵}_B$.

Finally, Figure 9 shows how the power generated by the different controllers varies according to the uncertainty added to the hydrostatic restoring coefficient of the system. In this aspect, it can be appreciated how the robustness of all the controllers is more limited. If the value of the coefficient $k_{res}$ increases, the excitation force that the wave exerts on the system (53) increases proportionally. Therefore, the margin of action of the PTO system decreases noticeably. Even so, the MPC₅ supports an added uncertainty of $25\%$, again being the one that provides the best robustness results even with the smallest prediction time $t_f$. Note that, for such added uncertainty, the excitation force becomes more than three-times the force that the PTO system can apply to the buoy; see Figure 10. After the MPC₅, the MPC₄ and MPC₁ get the best results, in this order.

5. Conclusions

The interest in implementing MPCs in WEC systems is motivated by the need to increase the productive/economic viability of these systems. For this reason, in this work, five different predictive controllers have been designed. All these controllers allow minimizing mechanical fatigue by limiting the operating range of the WEC system by means of hard and soft constraints. The main contribution of this work is the study of performance and robustness carried out for the five MPCs designed. This study demonstrates how the addition of an embedded integrator to the design model improves the performance of the WEC device referring to: average power generated, diminution of instantaneous power peaks, reduction of mechanical fatigue and robustness of the closed-loop system, in comparison with the other MPCs. In addition, the MPCs have been compared with two conventional controllers for WEC systems (resistive damping and I-P control) in a realistic sea state defined by the JONSWAP spectrum. Predictive controllers have proven that, compared to these conventional controllers, they minimize mechanical fatigue due to overshoot of nominal limits and increase the mechanical power generated by this type of WEC system.