Next Article in Journal
A Real-Time Pricing Scheme for Energy Management in Integrated Energy Systems: A Stackelberg Game Approach
Next Article in Special Issue
Coupling Methodology for Studying the Far Field Effects of Wave Energy Converter Arrays over a Varying Bathymetry
Previous Article in Journal
Coupling Analysis and Performance Study of Commercial 18650 Lithium-Ion Batteries under Conditions of Temperature and Vibration
Previous Article in Special Issue
Assessing the Macro-Economic Benefit of Installing a Farm of Oscillating Water Columns in Scotland and Portugal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Escuela Superior de Ingeniería, Universidad de Cádiz, 11519 Puerto Real, Spain
*
Author to whom correspondence should be addressed.
Energies 2018, 11(10), 2857; https://doi.org/10.3390/en11102857
Submission received: 27 September 2018 / Revised: 18 October 2018 / Accepted: 19 October 2018 / Published: 22 October 2018
(This article belongs to the Special Issue Offshore Renewable Energy: Ocean Waves, Tides and Offshore Wind)

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 z ¨ = F e + F r + F s + F u
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 and F r K r . The first is a function of the acceleration of the system and the mass of water added m ; while the second defines a radiation force as a function of the speed of oscillation as an integral of convolution.
F r = m w ˙ ( t ) F r m t h r ( t τ ) w ( τ ) d τ F r K r
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).
F r K r ( s ) W ( s ) = K r ( s ) = a n s n + a n 1 s n 1 + + a 0 b n s n + b n 1 s n 1 + + b 0
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 r e s = ( V d e s p ( z ) ρ g m g )
F r e s = k r e s z ( t ) k r e s = ρ g A W
where V d e s p is the volume of water displaced, ρ is the density of seawater, g is the gravity constant, k r e s 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 F e ( s ) = F e ( s ) η ( s ) = K τ s 2 + 2 ζ τ ω n τ s + ω n τ 2 .
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 η ¨ ( t ) + B a p r o x η ˙ ( t ) + k r e s η ( t )
where η is the height of the incoming wave and B a p r o x 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 p t o ), is given by (8).
G p t o ( s ) = F p t o ( s ) F u ( s ) = d a ω n 2 s 2 + 2 ζ ω n s + ω n 2 ,

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 r e s . 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.
F e ( s ) η ( s ) = 6.5 × 10 5 ( s + 0.9 ) 2 N m
K r ( s ) = F r ( s ) W ( s ) = 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 N m / s
where the coefficients a K and b K are listed in Table 1.
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 W E C ( s ) = Z ( s ) F e x t ( s ) = 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 m N
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.
x ˙ r ( t ) = A r x r ( t ) + B r w ( t ) F r K r ( t ) = C r x r ( t ) + D r w ( t )
where the matrices A r , B r , C r and D r are given by (13).
A r = 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 , B r = 121.3816 1.3375 0.1201 0.0005 C r = 10 3 1.0083 0.3683 0.2624 0.0377 , D r = 1218.70
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 p t o = 8.7965 1.0000 157.9137 0 , B p t o = 0 157.9137 , C p t o = 1 0 , D p t o = 0
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 p t 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.
z ˙ w ˙ x ˙ r x ˙ p t o x ˙ = 0 1 0 0 k r e s m T 1 m T D r 1 m T C r 1 m T C p t o 0 B r A r 0 0 0 0 A p t o A W E C z w x r x p t o x + 0 0 0 B p t o B W E C u F u + 0 1 m T 0 0 B W E C F e F e z w y = 1 0 0 0 0 1 0 0 C W E C z w x r x p t o x
where m T = m + m 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.

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.
z ˙ w ˙ x ˙ = 0 1 k r e s m T B a p r o x m T A W E C r z w x + 0 1 m T B W E C r ( F u + F e ) z w y = 1 0 0 1 C W E C r z w x
where m T = m + m , 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. M P 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 = A A 2 A 3 A M J x ( M n × n ) x k + B 0 0 0 A B B 0 0 A 2 B A B B 0 A M 1 B A M 2 B A M 3 B A M N B J u ( M n × N ) ( F p t o + F e )
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 p t o 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 p t o + F e )
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 g e n ( t ) = w ( t ) F p t o ( t )
P g e n = W T F p t o
where W and F p t o are vectors with length M, which represent the oscillation speed and control force for the whole prediction horizon, respectively.
By replacing (18) in (20),
P g e n = ( S w X ) T F p t o P g e n = ( S w ( J x x k + J u F p t o + J u F e ) ) T F p t o
where S w is a selector matrix for speed w (size M × M n ).
Developing (21) and grouping in terms of least squares, the cost function is obtained:
J ( F p t o ) = 1 2 F p t o T ( J u T S w T + R ) H F p t o + 1 2 ( F e T J u T S w T + x k T J x T S w T ) b F p t o
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:
J ( F p t o ) = 1 2 F p t o T H F p t o + 1 2 b F p t o A g F p t o B g
where A g = [ A 1 A 2 A 3 ] T and B g = [ B 1 B 2 B 3 ] T ; see Equation (24).
I I A 1 F p t o F P T O m a x F P T O m i n B 1 , S z J u S z J u A 2 F p t o Z n m a x S z ( J x x k + J u F e ) Z n m i n + S z ( J x x k + J u F e ) B 2 S w J u S w J u A 3 F p t o W n m a x S w ( J x x k + J u F e ) W n m i n + S w ( J x x k + J u F e ) B 3
where S z is a selector matrix for position (size M × M n ), the vectors F P T O m a x and F P T O m i n (size N × 1 ) define the nominal limits of the force applied by the PTO and the vectors W n m a x , W n m i n , Z n m a x and Z n m i n (size M × 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 p t o , ε z , ε w ) = 1 2 F p t o T H F p t o + b F p t o + ε z T W ε z ε z + ε w T W ε w ε w
where ε z and ε w represent the relaxation applied to position and speed along the prediction horizon M and the matrices W ε z and W ε w (size M × M ) weight these slacks.
By regrouping terms and adding soft restrictions, the cost function (25) can be expressed as:
J ( β ) = 1 2 β T H 0 0 0 W ε z 0 0 0 W ε w β + b 0 0 β A g β B g
where β = [ F p t o ε z ε w ] T , with size ( N + 2 M ) × 1 , A g = [ A 1 A 2 A 3 A 4 A 5 ] T and B g = [ B 1 B 2 B 3 B 4 B 5 ] T ; see Equation (27). Matrices 0 have the required size according to their location.
S z J u I 0 S z J u I 0 A 2 β Z n m a x S z ( J x x k + J u F e ) Z n m i n + S z ( J x x k + J u F e ) B 2 , S w J u 0 I S w J u 0 I A 3 β W n m a x S w ( J x x k + J u F e ) W n m i n + S w ( J x x k + J u F e ) B 3 I 0 0 I 0 0 A 1 β F P T O m a x F P T O m i n B 1 , 0 I 0 0 I 0 A 4 β κ z 0 B 4 , 0 0 I 0 0 I A 5 β κ w 0 B 5
where κ z and κ w are vectors (size M × 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. M P 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 r e f ) along a prediction horizon M,
J = ( w ˜ w r e f ) T Q ( w ˜ w r e f ) + F p t o T R F p t o
where Q and R are diagonal matrices of size ( M × M ) and ( N × 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 p t o ) = ( J x x k + J u F p t o + J u F e f w r e f ) T Q ( J x x k + J u F p t o + J u F e f w r e f ) + F p t o T R F p t o
By developing the cost function (29) and grouping terms, the following expression is obtained:
J ( F p t o ) = F p t o T ( J u T δ J u + R ) H F p t o + 2 ( f w r e f ) T Q J u b F p t o + ( f w r e f ) T Q ( f w r e f ) l
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 p t o ),
J ( F p t o ) = 1 2 F p t o T ( J u T δ J u + R ) H F p t o + ( J x x k + J u F p t o + J u F e w r e f ) T Q J u b F p t o
The constraints imposed on this cost function can be expressed in the same way as in the M P 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 r e f = F e / 2 B a p r o x .

3.3. M P 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:
x ( t + 1 ) y ( t + 1 ) x e ( t + 1 ) = A 0 C A I A e x ( t ) y ( t ) x e ( t ) + B C B B e F p t o ( t ) + F e ( t ) y ( t ) = 0 I C e x ( t ) y ( t ) x e ( t )
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 = C A C A 2 C A 3 C A M F ( 2 M × n e ) X e + C B 0 0 0 C A B C B 0 0 C A 2 B C A B C B 0 C A M 1 B C A M 2 B C A M 3 B C A M N B G ( 2 M × N ) ( F p t o + F e )
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 p t o ) for the full control horizon N,
J = ( w ˜ w r e f ) T Q ( w ˜ w r e f ) + F p t o T R F p t o
where Q and R are diagonal matrices of size ( M × M ) and ( N × 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 p t o :
J ( F p t o ) = 1 2 F p t o T ( G T δ G + R ) H F p t o + ( F x e k + G F e w r e f ) T Q G b F p t o
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:
J ( F p t o ) = 1 2 F p t o T H F p t o + b F p t o A g F p t o B g
where A g = [ A 1 A 2 A 3 ] T and B g = [ B 1 B 2 B 3 ] T ; see Equation (37).
T T A 1 F p t o F P T O m a x F P T O m i n B 1 , S z G S z G A 2 F p t o Z n m a x S z ( F x e k + G F e ) Z n m i n + S z ( F x e k + G F e ) B 2 S w G S w G A 3 F p t o W n m a x S w ( F x e k + G F e ) W n m i n + S w ( F x e k + G F e ) B 3
where S z and S w are selector matrices for z and w (size M × M n ), T is a lower triangular matrix (size M × N ), the vectors F P T O m a x and F P T O m i n (size N × 1 ) define the nominal limits of the force applied by the PTO and the vectors W n m a x , W n m i n , Z n m a x and Z n m i n (size M × 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,
J ( β ) = 1 2 β T H 0 0 0 W ε z 0 0 0 W ε w β + b 0 0 β A g β B g
where A g = [ A 1 A 2 A 3 A 4 A 5 ] T and B g = [ B 1 B 2 B 3 B 4 B 5 ] T ; see Equation (38).
T 0 0 T 0 0 A 1 β F P T O m a x F P T O m i n B 1 S z G I 0 S z G I 0 A 2 β Z n m a x S z ( F x e k + G F e ) Z n m i n + S z ( F x e k + G F e ) B 2 , 0 I 0 0 I 0 A 4 β κ z 0 B 4 S w G 0 I S w G 0 I A 3 β W n m a x S w ( F x e k + G F e ) W n m i n + S w ( F x e k + G F e ) B 3 , 0 0 I 0 0 I A 5 β κ w 0 B 5
where κ z and κ w are vectors (size M × 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. M P 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 M P 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 p t o + J f F e
where J x is a matrix already defined in Equation (17), while the matrices J u and J f are given by:
J u = B u 0 0 0 A B u B u 0 0 A 2 B u A B u B u 0 A M 1 B u A M 2 B u A M 3 B u A M N B u J f = B F e 0 0 0 A B F e B F e 0 0 A 2 B F e A B F e B F e 0 A M 1 B F e A M 2 B F e A M 3 B F e A M N B F e
where A, B u and B F e are obtained by discretizing (ZOH) A W E C , B W E C u and B W E 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 M P C 1 .
J ( F p t o ) = 1 2 F p t o T ( J u T S w T + R ) H F p t o + 1 2 ( F e T J f T S w T + x k T J x T S w T ) b F p t o
In addition, the nominal constraints imposed on the system must be added, which are defined in the same way as in the controller M P 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:
J ( F p t o ) = 1 2 F p t o T H F p t o + 1 2 b F p t o A g F p t o B g
where A g = [ A 1 A 2 A 3 ] T and B g = [ B 1 B 2 B 3 ] 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 M P 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:
J ( β ) = 1 2 β T H 0 0 0 W ε z 0 0 0 W ε w β + b 0 0 β A g β B g
where A g = [ A 1 A 2 A 3 A 4 A 5 ] T and B g = [ B 1 B 2 B 3 B 4 B 5 ] T ; see Equation (27).

3.5. MPC5

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 r e f . As in the M P 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 p t o + F x e k + G f F e f
where F is a matrix already defined in Equation (33), and the matrices G u and G f are given by:
G u = C e B e u 0 0 0 C e A e B e u C e B e u 0 0 C e A e 2 B e u C e A e B e u C e B e u 0 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 C e A e M N B e u G f = C e B e F e 0 0 0 C e A e B e F e C e B e F e 0 0 C e A e 2 B e F e C e A e B e F e C e B e F e 0 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 C e A e M N B e F e
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 M P C 3 , the cost function to be minimized can be expressed according to:
J ( F p t o ) = 1 2 F p t o T ( G u T δ G u + R ) H F p t o + ( F x e k + G f F e w r e f ) T Q G u b F p t o
In addition, it is necessary to add nominal constraints to the system, which are defined as in the controller M P 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:
J ( F p t o ) = 1 2 F p t o T H F p t o + 1 2 b F p t o A g F p t o B g
where A g = [ A 1 A 2 A 3 ] T and B g = [ B 1 B 2 B 3 ] 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 M P 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:
J ( β ) = 1 2 β T H 0 0 0 W ε z 0 0 0 W ε w β + b 0 0 β A g β B g
where A g = [ A 1 A 2 A 3 A 4 A 5 ] T and B g = [ B 1 B 2 B 3 B 4 B 5 ] 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 p t o ( k ) = K R D w ( k )
where k is the current sampling time and K R D 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 .
F p t o ( k ) = u P ( k ) + u I ( k ) u P ( k ) = K P w ( k ) u I ( k ) = K I T m ( w r e f ( k ) w ( k ) ) + u I ( k 1 )
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 p t o 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.

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.
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 (MPC3 and MPC5) achieve more regular power than the MPCs most used for WEC systems, those whose optimization criteria maximize the extracted power directly (MPC1 and MPC4). 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 MPC2 to the system, it gives the most irregular power; while the power generated with the MPC3, 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 MPC1 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 shows how the addition of the embedded integrator (MPC3) considerably improves the behavior of the system with respect to that obtained by applying the MPC2. Since, the MPC3 does not exceed nominal limits of the position and oscillation speed on any occasion. Furthermore, the MPC3 controller generates a clearer control signal than the MPC2 (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 MPC5 controller is the one that generates more power, followed very closely by the MPC1, MPC4 and, with a bit more distance, the MPC3. On the other side, the indicators O N L P and O N L S quantify the area of overshoot from the nominal limits of the position and oscillation speed, respectively. In this sense, the MPC3 (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 MPC1. 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 r e f 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 MPC1 and MPC4, 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 MPC3 and MPC1. 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 MPC3 does not have a monotonously increasing behavior that relates t f to the power generated, as would be expected. Nevertheless, in this aspect, the MPC3 is better than the MPC1, 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 MPC3 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 r e s (a nonlinear parameter that has been linearized during the modeling of the WEC system).
F r e s ( t ) = k r e s ( 1 + k r e s ) z ( t ) , F r m ( t ) = m ( 1 + m ) w ˙ ( t ) , F r K r ( s ) W ( s ) = ( 1 + B ) K r ( s )
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 ( 1 + m ) η ¨ ( t ) + B ( 1 + B ) η ˙ ( t ) + k r e s ( 1 + k r e s ) η ( t )
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.
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.
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 , the MPC2 is the least robust. Meanwhile, the MPC1 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 MPC5 extracts significantly more power than the others. It should also be noted that the controllers that directly maximize power in their cost function (MPC1 and MPC4) have the most predictable behavior with respect to the added uncertainty in m . On the other hand, the MPC5 offers the best features with respect to the uncertainty added to the dynamics of the radiation force (up to 400 % ). In contrast, the MPC2 gives very bad results in this respect. The MPC1 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 r e s 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 MPC5 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 MPC5, the MPC4 and MPC1 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.

Author Contributions

R.G. was responsible for designing the predictive controllers, studying them and writing most of the paper. A.C. and M.J.L. performed the modeling of the WEC system, proposed the identification methodology for PAW systems and the comparative robustness, reviewed all the work done and wrote part of the document. Furthermore, M.J.L. and A.C. proposed the control problems to solve using the model predictive control strategy and participated in the design and analysis of the controllers.

Funding

This project is funded by the Spanish Ministry of Economy, Industry and Competitiveness within the framework of the State Plan for Scientific and Technical Research and Innovation 2013-2016/State Programme for R+D+I aimed at the challenges of society.

Acknowledgments

Thanks to all ORPHEO (Optimización de la Rentabilidad de Plataformas Híbridas de energía Eólica y de las Olas) project partners for giving us the opportunity and funds to research in this promising field.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Pecher, A.; Kofoed, J.P. Handbook of Ocean Wave Energy; Sprinder Open: Boca Raton, FL, USA, 2016; ISBN 978-3-319-39888-4. [Google Scholar]
  2. The Ocean Energy Systems Technology Collaboration Programme. Available online: www.ocean-energy-systems.org/index.php (accessed on 4 May 2017).
  3. Pelagic Power. Available online: www.pelagicpower.no/about.html (accessed on 6 May 2017).
  4. Falnes, J. Ocean Waves and Oscillating Systems; Cambridge University Press: New York, NY, USA, 2002; ISBN 0-521-78211-2. [Google Scholar]
  5. Drew, B.; Plummer, A.R.; Sahinkaya, M.N. A review of wave energy converter technology. J. Power Energy 2009, 223, 887–902. [Google Scholar] [CrossRef] [Green Version]
  6. Valério, D.; Beirao, P.; Mendes, M.J.G.C.; da Costa, J.S. Comparison of control strategies performance for a Wave Energy Converter. In Proceedings of the 16th Mediterranean Conference on Control and Automation, Ajaccio, France, 25–27 June 2008; pp. 773–778. [Google Scholar]
  7. Li, G. Predictive control of a wave energy converter with wave prediction using differential flatness. In Proceedings of the 54th IEEE Conference on Decision and Control, Osaka, Japan, 15–18 December 2015; pp. 3230–3235. [Google Scholar]
  8. Andersen, P.; Pedersen, T.S.; Nielsen, K.M.; Vidal, E. Model Predictive Control of a Wave Energy Converter. In Proceedings of the 2015 IEEE Conference on Control Applications (CCA 2015), Sydney, Australia, 21–23 September 2015; pp. 1540–1545. [Google Scholar]
  9. Brekken, T.K. On Model Predictive Control for a point absorber Wave Energy Converter. In Proceedings of the IEEE Trondheim PowerTech, Trondheim, Norway, 19–23 June 2011; pp. 1–8. [Google Scholar]
  10. Richter, M.; Magana, M.E.; Sawodny, O.; Brekken, T.K.A. Nonlinear Model Predictive Control of a Point Absorber Wave Energy Converter. IEEE Trans. Sustain. Energy 2013, 4, 118–126. [Google Scholar] [CrossRef]
  11. Li, G.; Belmont, M.R. Model predictive control of sea wave energy converters—Part I: A convex approach for the case of a single device. Renew. Energy 2014, 69, 453–463. [Google Scholar] [CrossRef]
  12. Cavaglieri, D.; Bewley, T.R.; Previsic, M. Model Predictive Control leveraging Ensemble Kalman forecasting for optimal power take-off in wave energy conversion systems. In Proceedings of the American Control Conference, Chicago, IL, USA, 1–3 July 2015; Volume 2015, pp. 5224–5230. [Google Scholar]
  13. Oetinger, D.; Magaña, M.E.; Member, S.; Sawodny, O. Decentralized Model Predictive Control for Wave Energy Converter Arrays. IEEE Trans. Sustain. Energy 2014, 5, 1099–1107. [Google Scholar] [CrossRef]
  14. Li, G.; Belmont, M.R. Model predictive control of a sea wave energy converter: A convex approach. In Proceedings of the International Federation of Automatic Control (IFAC 2014), Cape Town, South Africa, 24–29 August 2014; Volume 19, pp. 11987–11992. [Google Scholar]
  15. Soltani, M.N.; Sichani, M.T.; Mirzaei, M. Model Predictive Control of Buoy Type Wave Energy Converter. In Proceedings of the International Federation of Automatic Control (IFAC 2014), Cape Town, South Africa, 24–29 August 2014; Volume 47, pp. 11159–11164. [Google Scholar]
  16. Starrett, M.; So, R.; Brekken, T.K.A.; McCall, A. Increasing power capture from multibody wave energy conversion systems using model predictive control. In Proceedings of the 2015 IEEE Conference on Technologies for Sustainability (SusTech), Ogden, UT, USA, 30 July–1 August 2015; pp. 20–26. [Google Scholar]
  17. Lagoun, M.S.; Benalia, A.; Benbouzid, M.E.H. A predictive power control of Doubly fed induction generator for wave energy converter in irregular waves. In Proceedings of the 1st International Conference on Green Energy, Sfax, Tunisia, 25–27 March 2014; pp. 26–31. [Google Scholar]
  18. Camacho, E.F.; Bordons, C. Model Predictive Control; Sprinder: Sevilla, Spain, 1998; ISBN 978-1-85233-694-3. [Google Scholar]
  19. Wang, L. Model Predictive Control System Desing and Implementation Using MATLAB; Sprinder: Melbourne, Australia, 2009; ISBN 978-1-84882-330-3. [Google Scholar]
  20. Openore. Available online: https://openore.org/2016/04/28/openwec/ (accessed on 14 May 2017).
  21. Ricci, P. Time-Domain Models. In Numerical Modelling of Wave Energy Converters; Elsevier Inc.: New York, NY, USA, 2016; pp. 31–66. ISBN 978-0-12-803210-7. [Google Scholar]
  22. Bozzi, S.; Miquel, A.M.; Antonini, A.; Passoni, G.; Archetti, R. Modeling of a point absorber for energy conversion in Italian seas. Energies 2013, 6, 3033–3051. [Google Scholar] [CrossRef] [Green Version]
  23. Hong, Y.; Eriksson, M.; Boström, C.; Waters, R. Impact of generator stroke length on energy production for a direct drivewave energy converter. Energies 2016, 9, 730. [Google Scholar] [CrossRef]
  24. Hansen, R.H.; Kramer, M.M. Modelling and Control of the Wavestar Prototype. In Proceedings of the 9th European Wave and Tidal Energy Conference, Southampton, UK, 5–9 September 2011; pp. 1–10. [Google Scholar]
  25. Kovaltchouk, T.; Multon, B.; BenAhmed, H.; Glumineau, A.; Aubry, J. Influence of control strategy on the global efficiency of a Direct Wave Energy Converter with electric Power Take-Off. In Proceedings of the Eighth International Conference and Exhibition on Ecological Vehicles and Renewable Energies, Monte Carlo, Monaco, 27–30 March 2013; pp. 1–10. [Google Scholar]
  26. Fusco, F.; Ringwood, J.V. A study of the prediction requirements in real-time control of wave energy converters. IEEE Trans. Sustain. Energy 2012, 3, 176–184. [Google Scholar] [CrossRef]
  27. Tedeschi, E.; Carraro, M.; Molinas, M.; Mattavelli, P. Effect of control strategies and power take-off efficiency on the power capture from sea waves. IEEE Trans. Energy Convers. 2011, 26, 1088–1098. [Google Scholar] [CrossRef]
  28. Cummins, W. The Impulse Response Function and Ship Motions. Schiffstechnick 1962, 9, 101–109. [Google Scholar]
  29. Morison, J.R.; O’Brien, M.P.; Johnson, J.W.; Schaaf, S.A. The forces exerted by surface waves on piles. Soc. Pet. Eng. 1950, 189, 149–154. [Google Scholar] [CrossRef]
  30. Openore. Available online: https://openore.org/2014/01/21/nemoh-open-source-bem/ (accessed on 13 October 2018).
  31. Armesto, J.A.; Guanche, R.; Jesus, F.; Iturrioz, A.; Losada, I.J. Comparative analysis of the methods to compute the radiation term in Cummins’ equation. J. Ocean Eng. Mar. Energy 2015, 1, 377–393. [Google Scholar] [CrossRef] [Green Version]
  32. Duarte, T.; Sarmento, A.; Alves, M.; Jonkman, J. State-Space Realization of the Wave-Radiation Force within FAST. In Proceedings of the ASME 2013 32nd International Conference on Ocean, Offshore and Arctic Engineering, Nantes, France, 9–14 June 2013; pp. 1–12. [Google Scholar]
  33. Mathworks. Available online: https://es.mathworks.com/help/signal/ref/prony.html (accessed on 2 February 2017).
  34. Ogata, K. Ingeniería de Control Moderna; Pearson Educación S.A.: Madrid, Spain, 2010; ISBN 9788483226605. [Google Scholar]
  35. Mathworks. Available online: https://es.mathworks.com/help/optim/ug/quadprog.html (accessed on 17 May 2017).
  36. Fusco, F.; Ringwood, J.V. Short-Term Wave Forecasting for Real-Time Control of Wave Energy Converters. IEEE Trans. Sustain. Energy 2010, 1, 99–106. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Illustrative scheme: point absorber wave (PAW) system of the W2POWER platform.
Figure 1. Illustrative scheme: point absorber wave (PAW) system of the W2POWER platform.
Energies 11 02857 g001
Figure 2. (a) Mesh defined by openWEC for a cylindrical buoy (11 m long and 7.5 m in diameter). Comparison of the identified transfer functions: (b) Bode diagram for excitation force F e ( s ) ; (c) impulse response for radiation force K r ( s ) .
Figure 2. (a) Mesh defined by openWEC for a cylindrical buoy (11 m long and 7.5 m in diameter). Comparison of the identified transfer functions: (b) Bode diagram for excitation force F e ( s ) ; (c) impulse response for radiation force K r ( s ) .
Energies 11 02857 g002
Figure 3. Comparison of models: complete ( z 1 , w 1 ) vs. simplified ( z 2 , w 2 ). The height of the wave is n.
Figure 3. Comparison of models: complete ( z 1 , w 1 ) vs. simplified ( z 2 , w 2 ). The height of the wave is n.
Energies 11 02857 g003
Figure 4. Simulation of the instantaneous mechanical powers generated by the WEC system when applying the seven controllers to the mathematical model (15).
Figure 4. Simulation of the instantaneous mechanical powers generated by the WEC system when applying the seven controllers to the mathematical model (15).
Energies 11 02857 g004
Figure 5. Simulation of: wave force, force setpoint demanded for the PTO and real force produced by the PTO, when applying the seven controllers to the mathematical model (15).
Figure 5. Simulation of: wave force, force setpoint demanded for the PTO and real force produced by the PTO, when applying the seven controllers to the mathematical model (15).
Energies 11 02857 g005
Figure 6. Simulation of the positions obtained in the WEC system when applying the seven controllers to the mathematical model (15).
Figure 6. Simulation of the positions obtained in the WEC system when applying the seven controllers to the mathematical model (15).
Energies 11 02857 g006
Figure 7. Simulation of the oscillation speeds obtained in the WEC system when applying the seven controllers to the mathematical model (15).
Figure 7. Simulation of the oscillation speeds obtained in the WEC system when applying the seven controllers to the mathematical model (15).
Energies 11 02857 g007
Figure 8. Variation of the extracted mechanical power as a function of the prediction time for M P C 1 and M P C 3 (control parameters listed in Table 4).
Figure 8. Variation of the extracted mechanical power as a function of the prediction time for M P C 1 and M P C 3 (control parameters listed in Table 4).
Energies 11 02857 g008
Figure 9. Variations of the average mechanical power generated as a function of the added uncertainty to: added mass, damping coefficient and hydrostatic restoring coefficient.
Figure 9. Variations of the average mechanical power generated as a function of the added uncertainty to: added mass, damping coefficient and hydrostatic restoring coefficient.
Energies 11 02857 g009
Figure 10. Wave force, force setpoint demanded and real force produced by the PTO, by applying the MPC5 to the model (15), which has an added uncertainty to the hydrostatic restoring coefficient of 25 % .
Figure 10. Wave force, force setpoint demanded and real force produced by the PTO, by applying the MPC5 to the model (15), which has an added uncertainty to the hydrostatic restoring coefficient of 25 % .
Energies 11 02857 g010
Table 1. Model coefficients identified for the WEC system.
Table 1. Model coefficients identified for the WEC system.
CoefficientValueCoefficientValueCoefficientValueCoefficientValue
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
Table 2. Parameters of the models (15) and (16) for the WEC system.
Table 2. Parameters of the models (15) and (16) for the WEC system.
SymbolDescriptionValue
k r e s Hydrostatic restoring coefficient809,325 N/m
B a p r o x Stationary approximation to system damping21,497 N s/m
mMass of water displaced by the buoy at rest241,601.9 kg
m Mass of water added167,700 kg
Table 3. Hard and soft constraints for the WEC system.
Table 3. Hard and soft constraints for the WEC system.
SymbolDescriptionValue
F P T O m a x Maximum stationary force for the power take-off system450 K N
F P T O m i n Minimum stationary force for the power take-off system−450 K N
z n m a x Maximum nominal limit for the buoy position1.25 m
z n m i n Minimum nominal limit for the buoy position−1.25 m
w n m a x Maximum nominal limit for oscillation speed1 m/s
w n m i n Minimum nominal limit for oscillation speed−1 m/s
z f m a x Maximum physical limit for the buoy position1.7 m
z f m i n Minimum physical limit for the buoy position−1.7 m
w f m a x Maximum physical limit for oscillation speed1.3 m/s
w f m i n Minimum physical limit for oscillation speed−1.3 m/s
κ z Maximum slack applied to the position nominal limit0.45 m
κ w Maximum slack applied to the speed nominal limit0.3 m/s
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).
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).
ControllerRQ W ε z W ε w K P K I
MPC14.500 × 10−71.000 × 10101.000 × 107
MPC21.000 × 10−79.150 × 1041.000 × 1072.000 × 109
MPC35.000 × 1054.575 × 1045.000 × 10101.000 × 105
MPC44.500 × 10−71.000 × 1081.000 × 104
MPC55.000 × 10−54.000 × 1051.000 × 1071.000 × 109
RD7.902 × 105
I-P7.050 × 1052.228 × 104
Table 5. Results obtained in the application of the seven controllers to the WEC system. The powers are expressed in kW. Note that O N L P indicates overshoot of nominal limits for position and O N L S indicates overshoot of nominal limits for speed.
Table 5. Results obtained in the application of the seven controllers to the WEC system. The powers are expressed in kW. Note that O N L P indicates overshoot of nominal limits for position and O N L S indicates overshoot of nominal limits for speed.
Controller P ¯ gen ONLP ONLS P gen Max P gen Min
MPC1127.600.00000.0023437.48−356.33
MPC2110.720.00000.0281744.06−364.97
MPC3125.530.00000.0000444.24−183.62
MPC4127.500.00000.0107452.75−369.27
MPC5129.010.00000.0413481.53−182.85
RD67.020.00000.0000258.88−1.05
I-P109.900.00730.0511583.73−107.34
Table 6. Sinusoidal components used for sea-state (values expressed in international units).
Table 6. Sinusoidal components used for sea-state (values expressed in international units).
ComponentAmplitudePeriodPhaseComponentAmplitudePeriodPhase
s 1 0.420 13.00 π s 9 0.200 9.00 0.00
s 2 0.520 12.50 1.5 π s 10 0.180 8.50 π
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 π s 13 0.100 5.50 0.5 π
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
Table 7. Control parameters set for the sea state recorded in Table 7.
Table 7. Control parameters set for the sea state recorded in Table 7.
Controller T m [ s ] T f [ s ] RQ W ε z W ε w
MPC1 0.05 1.5 1.1 × 10 7 1.0 × 10 10 1.0 × 10 7
MPC2 0.04 1.2 1.0 × 10 7 2.15 × 10 4 1.0 × 10 7 2.0 × 10 9
MPC3 0.04 1.2 5.0 × 10 5 1.475 × 10 4 5.0 × 10 10 1.0 × 10 5
MPC4 0.05 1.5 4.5 × 10 7 1.0 × 10 8 1.0 × 10 4
MPC5 0.02 0.6 5.0 × 10 5 1.75 × 10 4 1.0 × 10 7 1.0 × 10 9

Share and Cite

MDPI and ACS Style

Guardeño, R.; Consegliere, A.; López, M.J. A Study about Performance and Robustness of Model Predictive Controllers in a WEC System. Energies 2018, 11, 2857. https://doi.org/10.3390/en11102857

AMA Style

Guardeño R, Consegliere A, López MJ. A Study about Performance and Robustness of Model Predictive Controllers in a WEC System. Energies. 2018; 11(10):2857. https://doi.org/10.3390/en11102857

Chicago/Turabian Style

Guardeño, Rafael, Agustín Consegliere, and Manuel J. López. 2018. "A Study about Performance and Robustness of Model Predictive Controllers in a WEC System" Energies 11, no. 10: 2857. https://doi.org/10.3390/en11102857

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop