Sea trial results of a predictive algorithm at the Mutriku Wave power plant and controllers assessment based on a detailed plant model

Improving the power production in wave energy plants is essential to lower the cost of energy production from this type of installations. Oscillating Water Column is among the most studied technologies to convert the wave energy into a useful electrical one. In this paper, three control algorithms are developed to control the biradial turbine installed in the Mutriku Wave Power Plant. The work presents a comparison of their main advantages and drawbacks first from numerical simulation results and then with practical implementation in the real plant, analysing both performance and power integration into the grid. The wave-to-wire model used to develop and assess the controllers is based on linear wave theory and adjusted with operational data measured at the plant. Three different controllers which use the generator torque as manipulated variable are considered. Two of them are adaptive controllers and the other one is a nonlinear Model Predictive Control (MPC) algorithm which uses information about the future waves to compute the control actions. The best adaptive controller and the predictive one are then tested experimentally in the real power plant of Mutriku, and the performance analysis is completed with operational results. A real time sensor installed in front of the plant gives information on the incoming waves used by the predictive algorithm. Operational data are collected during a two-week testing period, enabling a thorough comparison. An overall increase over 30% in the electrical power production is obtained with the predictive control law in comparison with the reference adaptive


Introduction
Oceans worldwide cover around 70% of our planet. They concentrate a massive amount of energy, waiting to be harnessed. The challenge is to convert this energy effectively and in a reliable manner in this harsh environment. The global wave energy resource accessible along the globe coastlines is estimated to be between 0.5 and 2.2 TW [1]. Aware of that opportunity, scientists have developed hundreds of concepts to convert potential energy from the waves into electricity [2,3]. Among these Wave Energy Converters (WEC), only a small number have reached sea trials and today a few prototypes are being tested in real environment conditions. The high costs related to the development of such technologies and the high perceived risk tend to slow down the emergence of a sustainable market for wave energy conversion. Among the many existing concepts, the Oscillating Water Column (OWC) benefits a significant interest studied in several technology reviews [4e6]. Few prototypes reached sea trials and are actually in operation, or have been in operation for several years. To cite only the latest developments: the Pico plant [7e10], the Mutriku plant [11,12], the backward bend duct buoy from Ocean Energy [13e15], the Oceantec/Idom Marmok-A5 spar-type buoy [16,17]. One of the key advantages of the OWC is its simplicity and robustness due to its intrinsic working principle. Many existing PTO concepts, such as hydraulic PTO system or direct drive mechanical PTO, convert the slow motion of waves into mechanical energy and produce huge torque magnitudes when absorbing wave energy. The mechanical structure and components needs to be thought out to withstand these forces and can result in costly designs. In the OWC, the air trapped in a chamber adds another conversion step between the potential energy from the wave to the mechanical one at the turbine. This pneumatic power induced by the compression and expansion of air in the chamber creates an air flow that drives an air turbine rotating at high rotational speed and relatively low torque. The torque produced by the so-called self-rectifying turbines is unidirectional regardless of the flow direction and thus off-theshelf generators are suited for this type of turbines. The energy conversion chain of the OWC is displayed in Fig. 1. The wave energy device used in this study is one chamber of the Mutriku Wave Power plant (MWPP), that is situated in the North shore of Spain in the Bay of Biscay. It is located onshore and integrated into a breakwater. The total plant capacity using the 16 air chambers can reach up to 296 kW. However the first and the last chambers are disabled and thus the plant installed rated power totalises 260 kW thanks to 14 operational chambers, each one equipped with a set of Wells turbine and generator with 18.5 kW of nominal power [12]. A picture of the plant is presented in Fig. 2 and a cross-section of the chamber equipped with the biradial turbine is in Fig. 3. In operation since 2011, the plant is connected to the grid where it delivers the produced power. In the framework of the H2020 European funded project OPERA, the Mutriku power plant is at the centre of research activities. In this scope advanced control strategies are tested on a real environment [18] and specifically for the biradial turbine [19].
Among all the research area designated to reduce the cost of electrical production from wave energy sources, the discipline of control theory is focused in optimising the energy production. In the meantime, it requires very little additional investments being process instrumentation and control units. The interest given to the OWC technology resulted in a number of publications oriented to its control and gathered in a number of review publications [20,21]. Mathematical development concerning a fixed OWC was made for the first time by Nichols and Crossley [22]. In this research, optimal control techniques based on the projected gradient method were used for several strategies (pressure difference, flow rate and relief valve) to maximise the energy produced by the turbine. Numerical results were obtained with a model approximating the Pico plant applying optimal control to a simplified Power Take-Off (PTO) model. Improvements were seen to be within few percent to 35%. It arose phenomena to be taken into account but leading to unrealistic applications: chattering in valve operation, turbine allowed to be reversed to act as a compressor, the constrained pressure control was getting closer to the uncontrolled case while reducing the critical pressure. An air flow control was then implemented by Falcao and Justino [23] in a more complex numerical model of the plant including air compressibility, equivalent to a spring-like effect of the air chamber. Both the by-pass valve and series throttle valve were investigated and the effect on the production analysed for various turbine dampings. The benefits were quite substantial especially in avoiding the stalling behaviour of the Wells turbine. In Ref. [24] a stochastic modelling approach for the optimisation of the plant performance of the fixed OWC at the Pico Island equipped with a Wells turbine is proposed for a variety of realistic Sea States (SS). The aim was to optimise the turbine diameter and its optimal constant rotational speed that maximised the annual plant average efficiency with and without air flow control. The best improvement in terms of performance was found to be when optimising the turbine speed over using the flow control. Although a numerous examples of phase control strategies by using latching control have been studied [25e28], the present study focuses on turbine speed control strategies. In this line, numerous publications are focusing this type of control strategies. The FP7 european project CORES [14,29] leaded to experimental development of speed control algorithms for a Wells turbine and an impulse turbine with movable guide vanes to be tested in the 1/4th scale prototype of the OE buoy [30e32]. In total five algorithms were developed and implemented in a Hardware-In-the-Loop (HIL) electrical infrastructure. They were based on either pressure measurements or turbine speed to   compute the control torque to be applied by the generator. Also Henriques [33] developed several torque laws to control the biradial turbine installed in an OWC sparbuoy including a peak-power control that avoided generator overloads [34]. Again the turbine speed was controlled by applying an electromagnetic torque to the generator and these strategies were validated by HIL experiments. In Refs. [8,9] the Pico plant was used as a laboratory to develop and implement several algorithms to control the Wells turbine linked to the 200 kW generator. A flow control was permitted by a relief valve to avoid stalling issues produced by the turbine for high flow rates. An Artificial Neural Network algorithm estimated future potential turbine stall events enabling the valve operation to prevent the turbine from stalling. Several control strategies were tested on site and the performance analysis stated that the generator was oversized and the Wells turbine poorly characterised. It was found that besides the many assumptions used in the numerical simulations, the turbine characteristic curves were unrealistic. Other publications cites the Mutriku Wave power plant installation for the numerical development of control algorithms [35e39] to control the Wells turbine coupled to a Doubly-Fed Induction Generator (DFIG). Although they permitted to avoid stall behaviour of the turbine and improve its efficiency the model did not include the hydrodynamics of the plant meaning that the effect of the control on the internal pressure is not explored. An attempt to develop a more accurate model of the MWPP based on operational values coupled with wave data resource is presented in Ref. [40]. Although it is said that the 20 min average pressures are similar between the real and theoretical pressure drop, the time series of the pressure drop do not coincide in instantaneous values because of the smoothing effect of the approach. The model needed improvement to better represent the fast changes and all the dynamics of the real plant that is a very complex system. Lekube applied a Maximum Power Point Tracking controller to the model with the aim at controlling the Wells turbine speed and another flow control using the butterfly valve to avoid stall events [41,42]. The studies show an improvement of the power output with respect to an uncontrolled plant. However, the action of the change in the turbine speed on the internal water free surface was not considered. A hydrodynamic model of the plant including all the energy conversion steps was developed in Ref. [43] and used in Ref. [44] where Henriques proposed a comparison performance of the installed Wells turbine with the biradial turbine based on the common approach to design the control strategies.
The baseline of the present study relies in the Wave-to-Wire (W2W) model of the MWPP developed during the OPERA project [43]. The first contribution relies in using operational measurements and real wave climate data collected at the test site to improve the W2W model and make it more realistic. It serves as the common base to develop different control laws (CL) applied to the biradial turbine with its electrical conversion system. Thus, another contribution is the comparison of the CL, supported by numerical simulations, in terms of power production, power quality and practical implementation issues, taking into account the limits of the equipment used, i.e. the turbine and generator characteristics. All three control laws are variable speed control defined by their torque laws that are characterised using various approaches. One especially is based on a predictive algorithm made possible by the installation of a wave elevation sensor installed in front of the plant and communicating by wire with the instrumentation and control system onshore. The results show an improvement of the overall performance of the predictive algorithm as energy production is higher while the turbine and the generator operate within a safe range. The best adaptive controller and the model predictive control are then implemented in real environment in the wave power plant and evidence of operational results during field testing are presented to justify the improvements of the predictive controller in terms of power production. The datasets of results used to produce this paper are publicly available in the Zenodo platform in Ref. [45].
The paper is structured in seven sections organised as follows. Section 2 details the development of the numerical model including the methodology used to obtain accurate hydrodynamics frequency coefficients with operational data and create the W2W model. The air chamber model was also tuned thanks to data collected during the sea trials. Section 3 explains the power take-off system composed of the biradial turbine and a 30 kW induction generator. Then, section 4 develops the control strategies methodology and their implementation challenges. Section 5 deals with the numerical simulation results where the performance of each controller are assessed. Section 6 presents sea trial results and the improvements brought by the predictive algorithm are highlighted. Finally, the paper concludes the work in Section 7.

W2W model of one chamber of the MWPP
This section develops the methodology to better define the W2W numerical model. Data from measurements at the plant during operation are used to tune the model and improve its accuracy. More specifically the hydrodynamic coefficients are recalculated, and the air chamber model is adapted with a more realistic phenomenon.

Frequency domain model
The frequency domain analysis is done using a boundary element method solver based on potential flow theory and applied to the geometry of the OWC chamber of the real plant associated with the seabed slope in front the breakwater wall. Following the approach of the rigid piston modelling of the onshore OWC at Pico [7], the Internal Water Free Surface (IWFS) is considered a massless oscillating body restricted to the heave motion and the chamber wall is the second fixed body.
During the hydrodynamic characterisation it was found that both the seabed geometry and the length of the dike considered in the model would significantly affect the results obtained. This was already observed in Ref. [7], where the effects in the plant capture width was analysed. For the MWPP, and without detailed bathymetric information, the seabed geometry is simplified to a constant slope starting 50 m ahead of the plant with a 15 m water depth. The total length of the dike is 100 m (c.f. Fig. 4), where air chamber 9 would fall right in the centre, which results in a quite realistic assumption given the real plant layout. The hydrodynamic coefficients of excitation force F exc ðuÞ, radiation damping F rad ðuÞ and added mass at infinite frequency AðuÞ for each u angular frequency obtained from the hydrodynamic solver are represented in Fig. 5.
The excitation force appears higher in the low frequency range between 0.2 and 0.6 rad/s, where a maximum in force coefficients is observed around u ¼ 0:4 rad/s. In addition, the inclusion of a long dike (100 m) in the model introduces high frequency ripples in the excitation force, that otherwise don't exist with a short dike. The short dike would only consider the length of the two surrounding chambers at each side of the studied chamber. Moreover, if a very short dike length (or no dike length at all) would have been considered, the effect of the slope is negated. The results obtained are checked both with the Haskind relations and the pressure integration method [46], and the results obtained are quite consistent.

Internal water free surface motion
The W2W model in time domain, t being time, resolves the differential equation of motion in heave z of the IWFS following the 2nd law of Newton: where A ∞ is the added mass at infinite frequency of the OWC, € zðtÞ the heave acceleration, F exc ðtÞ the excitation force from the waves, F rad ðtÞ is the radiation force, and finally F p ðtÞ is the force of the air pressure. The excitation force of the incoming waves on the IWFS is expressed as: that is the sum of N wave components each one described as the product of F exc ðuÞ the excitation force from the frequency domain analysis with the individual wave amplitude A w and phase 4, lastly ε is the random phase from [0:2p]. The radiation force is the damping component due to the motion of the free surface over still water outside the wall and is obtained by solving the convolution integral equation: K is the impulse function that is approximated using the Prony method [47] (c.f. Fig. 6) with the coefficients a and b representing the matrices of the complex radiation approximation coefficients: The convolution integral function can therefore be approximated by P N i¼1 I i and mathematically solved by the differential equation:  A 5th order Prony approximation was found to be accurate enough as seen in Fig. 7 and the computation is solved by a substate-space system integrated in the global system.
The hydrostatic force is composed of the water density r w , the gravitational constant g, and S IWFS the area of the IWFS.
The differential equation is solved in time domain using monochromatic waves for each angular frequency without the pressure forces F p and the Response Amplitude Operator (RAO) of the heave motion is compared against the frequency domain model in order to validate the time domain W2W model. In Fig. 7, the difference between the two models is showed. It proves the correct implementation of the time domain model as the error is relatively small and falls in the range [5% to À10%], the largest errors occurring in the high frequencies only. An ideal reflection coefficient (¼1) is used when modelling the walls of the chamber and assuming that waves are totally reflected. This fact leads to a RAO value starting in a value above 1 for very low angular frequencies. Further work on the hydrodynamic model should consider more realistic values, where some reflections losses are allowed.

Air chamber model
The force F p , is the resultant of the pressure forces on the IWFS and is associated to the power extraction by the PTO. Considering p Ã the dimensionless relative pressure at the turbine inlet between p  the internal chamber pressure and p at the atmospheric pressure outside the chamber, the formulation is then p Ã ¼ p= p at À 1.
In order to define F p , it is convenient to develop an approach to model the air inside the chamber. Depending on the air chamber model, for example between one considering the air incompressible against another one where air is compressible, the magnitudes of pressure can be overestimated up to 15% [48] and have an impact on the output power between 10 and 30% [49]. In Ref. [50] the incompressible, unrealistic, air model overestimated the output power by a factor of 1.7 with respect to a model including air compressibility. The pressure variation has a non-linear spring-like effect on the internal water surface that has to be taken into account in the W2W model.
The model of air developed here considers air as a perfect gas subjected to a polytropic process. In Ref. [51], k is the polytropic exponent equivalent to the turbine polytropic efficiency, also known as small-stage efficiency, and considering g ¼ 1:4 the specific heat ratio: The turbine efficiency h t (c.f. Eq. (15)) is needed to compute the polytropic exponent. The pressure variation is expressed as: and depends on the turbine mass flow rate _ m t , the volume inside the chamber V ch ¼ S IWFS ðh 0 þ zÞ with respect to the air chamber height at rest h 0 . Here the air density of the chamber is set as the atmospheric one r ch ¼ r at . During a wave cycle, there is a cycle of compression and expansion of air inside the chamber. This implies a change in the air density r in at the turbine inlet and affects the turbine behaviour as explained in Ref. [50]. During compression the density is computed as: During the expansion phase, the air comes from outside the chamber at atmospheric conditions and is naturally r in ¼ r at . The variable air density behaviour is introduced in Section 3.1.
A comparison is done on two air compressibility models, one uses an isentropic process like the one in Ref. [43] and the other, a polytropic one. In fact, the isentropic process is a special case of the polytropic one when k ¼ g. When comparing the ranges of pressure the polytropic model shows lower amplitudes as illustrates Fig. 8, and seems more realistic. Eventually, the pressure variation differences have an impact on the power outputs given by Fig. 9 where the pneumatic, the turbine and the generator powers are displayed. Note that both figures refer to simulations results of SS10 (c.f. Fig. 18) where the control CL1 is employed (c.f. Section 5.2). The power quantities are explained in Section 3 and the first control method developed in Section 4.2 for this study. When comparing the mean power production between the polytropic model developed here and the isentropic one, it presents a decrease in electrical output power higher than 15%. In a variety of sea states, the difference of the three powers between the polytropic and the isentropic model can be seen in Fig. 10 and oscillates around À15%. The figure only presents the first 11 SS. Above, the sea states are too energetic and the safety valve begins to operate more often in the case using the isentropic model than in the other one.
The validation of the model based on the polytropic process is done with experimental data measured at the plant. The validation methodology consists in isolating the pressure variation. A numerical model takes as input the plant data of water surface elevation position and velocity, and the turbine rotational speed. The pressure variation obtained with the polytropic model is then compared with the measured pressure values. Fig. 11 shows the comparison of a modelled and measured relative dimensionless pressure for a medium energetic sea state. The model overestimates the pressure around 4% and the errors are more visible in frequencies lower to 0.8 rad/s. The isentropic model was not compared to the experimental data but the overestimation would have been higher. Following the same methodology for a range of realistic sea states where the comparison was undergone, the error is in the range of þ15% to À5%. Table 1 presents the main PTO characteristics for the biradial turbine and its generator installed in Mutriku. It is equipped with a High-Speed Shut-off Valve (HSSV) for safety use. The aerodynamic characteristics of the turbines used in the model can be represented by the relation F ¼ f ðJÞ, P ¼ f ðJÞ and h t ¼ f ðJÞ. The performance curves of the biradial turbine designed for OPERA are given in Fig. 12 [52].

Turbine model
In turbomachinery, it is common to define turbines according to the dimensionless aerodynamic parameters of pressure head J, mass flow rate F, power coefficient P and efficiency h t expressed as in Ref. [24]. As defined in Eq. (11), r in includes the change in air density during the compression phase and affects the turbine behaviour. Note the variable u v defines the safety valve position detailed in Section 4.1.

Generator model
The generator connected to the turbine shaft is a squirrel-cage induction one of 30 kW rated power and two pairs of poles, its main characteristics are gathered in Table

Loss model
The power conversion from mechanical to electrical is subjected to losses at the generator. The generator efficiency is defined by the amount of power extracted at the generator input and the sum of its power losses. The mechanical losses P ml are mainly due to bearing friction and depend on the diameter shaft and its number of poles, or nominal speed, and are estimated with data from the catalogue   [53] to be dozens of Watts. The core losses are related to eddy currents flowing in the laminated steel plates.
It depends on the magnetic field B, s h and s e are respectively the hysteresis and eddy current loss coefficients and their values are estimated through the generator presented in Ref. [54]. The steel plate thickness d ip is taken from Ref. [55]. The resistive losses from the winding circuit, mainly from copper in the stator, are calculated with the RMS stator current I rms and its resistance R st :   Finally the generator efficiency is so that h g ¼ P g =P in where the electrical power at the output of the generator is:

Operation ranges
The generator may operate for short periods of time over its rated capacity. The maximal torque that can be extracted at the nominal speed is set to be T max ¼ M mn T nom . This operation corresponds to the nominal voltage. When the turbine spins faster than the nominal, the voltage cannot be overshot and so the frequency raises above the nominal and there is still a torque extracted until the current limitation. This is the flux weakening region. This operation is represented in Fig. 13 but in the case studied in the paper the generator is connected to a power converter whose voltage rating is higher and its nominal capacity oversized. It means that the generator will operate normally until reaching the overspeed N os ¼ V pe =V g N nom ¼ 1:725 N nom . That arrangement allows to overload the generator close to three times the generator rated power. This is a challenging operation mode for the generator and thus requires a good ventilation in practice as the temperature inside the generator rises above normal condition. The model does take into account the maximal power extraction but does not include a heat model.

Drive train model
The angular acceleration, i.e. the sum of torques over the inertia of the turbo-generator set, is described such as: T t is the torque provided by the turbine and given by T t ¼ P t = U, T g is the resistive control torque to apply by the generator and T loss comes from the generator losses model.

Safety control of the PTO
A safety valve installed in front of the turbine will operate in high energetic sea states to protect the PTO components. It avoids over-speeding the turbine and reaching the maximum speed of the generator. If the threshold of the cut-off speed N co is reached, the valve instantaneously closes and blocks the air flow. This type of valve actually equips the biradial installed in Mutriku. In the closed position, the torque T g is applied at the generator following the control law and reduces the rotational speed until a cut-in speed N ci . When this value is reached, the valve opens and the turbine operates normally. The speed thresholds for the valve actuation on each turbine are given in Table 1

CL1: theoretical turbine torque estimation
A simple variable speed control strategy is designed to operate the turbine at its best operation point [31,56]. It is based on the assumptions that there is an optimal instantaneous torque provided by the turbine for a specific rotational speed (c.f. Best Efficiency Point (BEP) in Fig. 12. The objective is to control the turbine speed to reach this value by applying the resistive torque T g on the generator side. This control variable then affects the rotational acceleration in Eq. (19). For a lossless power conversion system the powers in the turbine and the generator side are equal in average. This means that, similarly to an optimal estimated turbine torque, there is an optimal torque at the generator side controller, defined by a quadratic law in function of the speed.
This kind of relation is called a torque-speed curve, or T-U curve, or torque law. The gain for the curve slope is selected for a specific design point offering the best turbine efficiency and is obtained rearranging equations (12)e(15): This T-U law presents the advantage to perform fair enough independently of the sea state and only needs turbine characteristics to compute the required design point. However it does not take into account the global dynamics of the plant nor the generator operation. It will serve as the base case to compare the simulation results of the controllers. The torque law is represented in Fig. 14. It appears that in the highest speeds, the controlled torque is only half of the nominal generator one.

CL2: optimum turbine operating region
This control strategy uses Eq. (20) to formulate it in a more general manner: It was employed in a number of publications [33,57,58]. The principle resides in finding torque law parameters best suited for a specific OWC operating in a certain wave climate and with a specific air turbine. It uses the numerical model to run an offline optimisation that finds the optimal operational condition. For each sea state, the methodology is to fix the turbine speed, by assuming an infinite inertia, upon Eq. (19) there is no acceleration, extract the average mechanical power produced by the turbine and repeat for a range of turbine speeds. For each sea state, there is an optimal speed associated to a maximum turbine average power. When performing a power law regression of the optimal operation points, the maximum achievable mechanical power is defined following Fig. 15. In the present case study, the torque law being associated to this operation conditions is defined by two control parameters being the T-U curve slope a ¼ 7e À 4 and its exponent b ¼ 2:1731. As for CL1, CL2 torque law is also plotted in Fig. 14. 4.4. CL3: non-linear constrained model predictive control

Predictive speed control
The aim of this non-linear model predictive control is to maximise the electrical output power. The algorithm is similar to the one presented in Ref. [56] the difference being that, in the present work, it is employed to a variable speed control instead of a pressure threshold latching strategy. The torque law proposed in the base case relies on a lossless drivetrain and on the perfect characterisation of the turbine. In comparison, the proposed strategy takes the overall power conversion by the PTO, and includes in the optimisation process the efficiencies of both the turbine and the generator. Eventually it includes all the numerical model non linearities. At each time step, the torque law in Eq. (22) is employed to compute the control torque, but using the manipulated variables of the curve slope a and its exponent b optimised by the controller. At the beginning of the control horizon, called replanning period, the best configuration of a and b is calculated for the prediction horizon T ph , and constitutes the new control vector u c ¼ ½a; b. These parameters are now updated and applied by the torque law at each time step during the replanning time, until a new control vector is computed at the beginning of the next replanning period. Fig. 16 depicts the timings of the controller.
The predictive horizon time is constrained by the wave travel time from the real-time wave elevation sensor to the plant.
Assuming shallow water conditions, the constant wave travel speed for all frequencies depends on the water depth d w , the theoretical wave celerity is defined as c w ¼ ffiffiffiffiffiffiffiffiffi gd w p . Because the sensor is installed at d s ¼ 200m ahead of the chamber wall at a mean tide depth around d w z10m, the average prediction time is This gives sufficient time to consider between one and two wave periods for the prediction. A replanning time of half the prediction time is appropriate in this case. As the optimisation takes into account average values of power, it is not so sensitive to synchronisation between the wave measured and the one actually hitting the plant. In fact in real conditions, the waves are not cylindrical and do not always travel in the same direction, adding uncertainty in the wave prediction.
The optimisation algorithm tends to minimise the cost function in Eq. (24). It includes the turbine and the generator power in order to integrate both efficiencies and constraints. Restriction on the generator maximal torque and the turbine threshold speed for the valve actuation are considered in the numerical model during the optimisation.
The parameters a J and b J are weighting parameters that allow setting priority to power transformation at the turbine or power conversion into electricity. A multi-parametric optimisation is made to find the best set of weights. The control parameters are able to evolve in respect to the previous values, by a maximum increment of ± 25%. This allows to guide the search of the optimal parameters during the optimisation.

Estimation of the predicted excitation force
One of the challenges of predictive algorithms is the estimation of the excitation force from the incoming waves. This issue has been raised in a number of publications [59e62] and the approach followed to resolve that issue is by applying the methodology of [63,64]. It consists in resolving the convolution between the Excitation Impulse Response Function (EIRF) K exc and the wave eleva- To do so the same method as the one used to approximate the radiation impulse function is applied but to the EIRF. The Prony approximation at order 5 gives fair results. The complex Prony matrices are integrated in a state-space system. The resolution of the estimated excitation force is done resolving this system having the wave elevation in input. The validation process in Fig. 17 is done by generating the excitation force and wave elevation respectively from Eqs. (2) and (26) with the one estimated using the present method. An error inside the range of ± 10% is obtained by comparing the two excitation forces modelled and estimated.

Simulation results
This section presents a quantitative comparison of the controllers performance. The numerical results are obtained running the W2W model with the three controllers and 14 sea states. The

The wave resource at Mutriku
The wave climate was measured by the pressure sensor Virtuoso from RBR, installed at the sea bottom in front of the plant during the period February to April 2018. Although there are uncertainties related to the measurements of the wave elevation using hydrostatic pressure variation, this equipment was the best trade-off between accuracy, reliability and cost. The wave elevation time series measured at 2 Hz were then used to obtain the 1 h sea states. These statistical wave data are available online on Zenodo [65], as for the methodology to extract the wave statistics and compute the significant wave height H s , and energy period T e . The selection of the 14 SS is done over the compromise of high occurrence for a large variety of sea states. Both the wave climate and the selected sea states are presented in Fig. 18. Each of the 14 wave spectra used in the simulations comes from all the real wave spectra inside the sea state bin and averaged over all the angular frequencies. Fig. 19 highlights the averaged spectra for all the SS. For example in Fig. 20 SS10 is the average of 90 real spectra measured in front of the plant.

Production performance
In this section the performance of the three control laws are compared. In Fig. 21, the equivalent annual energy production (AEP) of a CL is the sum of the energy produced during the 14 SS. It is computed by multiplying the average electrical power of a SS by the number of hours during which it operated, upon its occurrence during a year and assuming a 100% availability of the plant.
This figure gives a general view of the global power production performance. The improvement of CL2 in respect with CL1 is 23%. CL3 gives the highest score increasing the energy produced by nearly 50% from CL1 and 20% from CL2. Fig. 22 reveals the energy production divided by sea states. One can easily apprehend below a sea state of 1 m, from SS1 to SS6, there is not much difference in production between the control algorithms. From sea states higher, besides SS9 that rarely occurs, that is where the energy is produced. During mid-energetic sea states CL2 and CL3 are close and perform better than the base case. It is during more energetic sea states,  from SS10, that the predictive controller CL3 outruns the two other algorithms. Then, when comparing CL3 with CL2, the increase in energy production rises between 30% and 90%. During wave height above 2 m, the operation of the safety valve begins to have an effect on the power extraction at the turbine. CL1 and CL2 are not able to break the turbine enough to stay inside the speed operation range. Focusing on the electrical power production, Fig. 23 presents the performance of each CL per SS. Comparing CL2 and CL3, it appears that turnpoint begins at SS10 and the improvements of CL3 versus CL2 grow with the available wave power of the sea state. In the best cases, the first two CL are somehow clamped at an average power a third or half of the generator rated one. The explanation comes from the safety valve operation, the torque laws of CL1 and CL2 are too loose at high rotational speeds. CL3 in comparison is aware of the limitation on the turbine speed and can adapt the control parameters to produce steeper T-U curves. Fig. 24 shows the same figure as in Fig. 14 representing the CL torque laws but adding the CL3 simulation of SS10 associating the controlled torque to the rotational speed. This highlights the capacity of the predictive controller to adjust the torque law to the wave and plant conditions in a more dynamic manner. It is very versatile and can act as the CL1 for calm wave groups, and also apply a harder torque at high speeds to avoid turbine overspeeds and thus prevent the HSSV from staying closed during large time duration. Still, the generator limitations in terms of maximal torque are respected (c.f. Fig. 13). Table 2 allows a comparison of the valve actuation in the sea states where the HSSV operates by showing the time duration when the valve is closed. Obviously, the longer the valve stays closed, the higher the losses of absorbed power.
Analysing the results focusing on the efficiencies gives an additional perspective on the results analysis. In Fig. 25 the turbine efficiency is broken down for each SS. It reveals the average turbine efficiency is almost constant in all cases, removing the last SS, with    a little advantage of CL1 in the less energetic SS. Globally, CL3 appears here to perform slightly under the two others. The good score of CL1 is explained by the fact the controller operates the turbine at its highest efficiency point. A noticeable difference between the control strategies appears when analysing the generator efficiency.
Indeed, in all the sea states, CL3 is ahead of the others. This outcome is brought by considering the generator in the online optimisation. Finally, the total PTO efficiency, comprising the turbine and the generator ones, confirms the good score achieved by CL3 in terms of AEP. The improvements of the production when operating by CL3 are mainly due to better considerations given to the generator while the turbine operates almost the same for every CL.

Quality of power
This part highlights results in terms of power quality, in the sense of the amplitudes of power profile likely to create disturbances to the grid in larger production units or to affect the WEC electrical component life. Although at the scale of a 30 kW generator the impact of bad power quality to the grid is low, this study allows to understand how the controllers apprehend the peaks of production. The impacts and the stresses of power peaks on the electrical equipment caused by control strategies is out of the scope of this work but has been analysed in Ref. [66]. Mainly the standard deviation shows how far from the average electrical production are situated the extreme events. Also the peak-toaverage power ratio (Pk2av) is the ratio of the maximum value of electrical power from the generator over its average. Both of the two values quantifies the variation of power during a test. A high value represents a poorer quality of power injected in the grid. In Fig. 26, the Pk2av of all CL for each SS is represented. The figure reveals a much higher Pk2av for the predictive law against the others, sometimes doubling in absolute values. Globally, while the sea conditions get higher, the Pk2av reduces because in low sea states, the average electrical power is quite low. In order to complete the analysis, another way to value this quantity is proposed by Fig. 27. It shows the production peaks for each SS, normalised by the generator rated power, in function of the average generator load ratio L ¼ P g =P g;nom . CL1 presents peaks half of the rated power and mean powers, at most, reach a third of the rated power, a poor efficiency region for a generator. This control law exploits the generator way below its nominal capacity, the values of power peaks are harmless for the equipment. In CL2, the highest peaks correspond to the nominal generator power. Similarly to CL1, there is no generator overload. In the most powerful sea state, the generator power is at half its rated capacity. Same observation as the previous case, CL2 operates the generator in an inefficient manner. On the contrary, the predictive law is able to reach higher power peaks by applying steeper torque laws. In line with observing power peaks, the standard deviation of electrical power in Fig. 27 indicates the dispersion on the probability distribution. While CL1 and CL2 reach a stable deviation as the average load increases, CL3 follows a linear behaviour. The highest values of standard deviation reveal a larger variability of punctual high energy peaks coming from the fact that the predictive law is capable of absorbing more energetic group waves in high sea states. Consequently the power quality criteria drops for CL3. The power peaks reach twice or three times the rated generator power. That means in some cases the generator is operated at its limits. In the most energetic sea states, the average production rounds half the generator rated, and in the highest sea states, reaching almost the generator rated capacity. To contrast with the good score of power production, this law induces larger standard deviation and highest power peaks, that are most likely to reduce the generator life if it is often operated at these values. Induction generators can withstand, for short periods of time, high power peaks as long as the temperature of the windings does not overshoot a certain threshold. In a continuous overloaded operation the insulation of the copper windings can melt, cause a shortcut and seriously damage the generator. Although this operation represents a little percentage of the annual wave climate, in practice the control framework should include safe temperature thresholds that triggers alarms and adapts the power production in consequence. Also, either no results are presented, the predictive control can be easily tuned to reduce power peaks. This would consist in adding restrictions harder than the existing ones in the optimisation loop.
As a global analysis, at low to mid-energetic sea states CL2 and CL3 appear relatively close with 6% increase in energy production for CL3. Based on the simulation results, the advantage of CL3 are best perceived at energetic sea states starting from 2 m. In lower sea states it would be more questionable to use such a sophisticated controller instead of CL2. This controller may be preferred in areas where the resource is less important than in the location of Mutriku.

Sea trial results at Mutriku
During the OPERA project, the biradial turbine was installed in the 9th chamber of the MWPP during one year, from June 2017 to June 2018. Several control algorithms were tested and results are reported in the deliverable D4.2 [18]. During this testing period, the predictive controller CL3 could be tested and compared to a CL2 type of adaptive controller for a two week period at the end of June 2018. Unfortunately only low energetic sea states were experienced during the limited testing duration. A real time wave elevation sensor from Isurki was installed at the same point of the RBR Virtuoso (c.f. Section 5.1) and connected by cable to the plant instrumentation system to feed the predictive controller.

. Instrumentation and control framework
All the plant sensors were centralised and processed in the Programmable Logic Controller (PLC), a X20cCP1584 from B&R specially designed for harsh environmental conditions. Typical measurements included the motion of the IWFS, valves positions, pressures at different location of the turbine, the drive train speed, rotor axial vibrations and generator winding temperature, electrical quantities at the generator and the grid side of the Power Electronics (PE). All the operational data were sent by the PLC to a local database (DB) and replicated in a cloud web service for easy access and faster post processing. The global advisory and control software was designed following a state-machine approach with a security layer able to trigger different sorts of warnings or alarms. The hierarchy of alarms induced different safety actions to insure the plant and PTO integrity. The real-time control of the turbine speed was done by simply computing the torque reference and send it to the PE that follows this reference with its internal current and voltage control loops. The PLC cyclic time for the main control program was 100 ms. The test procedure was developed to assure autonomous operation and automatically switch the control laws between CL2 and CL3 each half an hour, the duration of a sea state.

Wave resource
The hydrostatic pressure sensor CNC4200-MT3 of Isurki was used for both providing the wave climate, as it was done with the RBR Virtuso, and sending the measures real time to the PLC. During the sea trials, the wave elevation sensor measured various SS but only 7 different sea states had statistical relevance considering a number of tests for both the CL. These SS are in the range of H s ¼ ½0:5 : 1 m and T e ¼ ½7 : 14 as seen in Fig. 28. These two weeks in June 2018 were poor in terms of variety of sea states but still 325 tests of half an hour each falling into that wave climate specification could be done. The real time data received by the PLC was sent as input to the predictive controller. The pre-processing included a tide cancellation to centre the wave elevation to 0 and estimated the excitation force. A moving vector, with a window of the prediction horizon time, was filled in and sent to the MPC at the beginning of each replanning period. To remove some uncertainty in the synchronisation between the MPC and the resource, the prediction time was variable and computed in function of the tide, thus the water depth, following a previous correlation study performed between the wave measured upwave by the RBR Virtuoso and plant operational data. The prediction time was between 20 and 30s, and as for the simulation, the replanning time was half the prediction one.

MPC implementation
To overcome the computation requirements of the on-line optimisation and the complexity in the implementation, it was decided to use an additional computer and link it to the PLC via an OPC server. This computer released the PLC from the computation burden and run a Matlab/Simulink model for the experiment of CL3. It received the wave elevation and the main plant operational data. At each time step, the excitation force moving vector was updated while the controller computed the radiation forces, as a virtual sensor would do, using the real IWFS velocity. At the beginning of each replanning, the MPC algorithm received the state vector initialised with the last instantaneous plant data and the computed radiation forces components, plus the vector of estimated prediction forces. The MPC computed the best pair of control parameters in a maximum of 3 PLC cyclic time and returned them to the PLC where they constituted the new torque law applied during a full replanning time. At the end of a 1/2 h test, the control parameters from CL3 were saved and initialised the next time the CL3 was called, after a 1/2 h operation of CL2, and so on.
During the implementation of the predictive algorithm, a 20-h test was done to validate its autonomous and safe operation. The evolution of control parameters in Fig. 29 is interesting in several ways. First it proves the controller was operational without interruption for a full day. In addition, at the beginning of this test, control parameters far from their optimal values were initialised to verify the algorithm could converge to a solution. Finally small variations of the a and b coefficients illustrates the fact the controller adapts to a short time changing wave conditions at each replanning period.

Power production performance
To facilitate the analysis and comparison of performance, the electrical power production is presented at first for CL2 over the 7 studied sea states in Fig. 30. Because the controllers were operational several times during a same sea state, a boxplot is a  convenient way to appreciate these test results. For each sea state, the mean of the averaged power production per sea state is shown with a star, the blue rectangle collects half of the tests between the quartiles 25% and 75%. The maximum and minimum are also represented at the extremities of the rectangles. The red crosses are outliers and thus not considered in the analysis. The dispersion of results is probably due to the fact that representing a sea state by H s and T e may not be the most appropriate for this case. Indeed, there are other criteria characterising a wave spectra such as the wave steepness or the peakedness of the spectrum. Also the wave direction and the tide level can influence the energy absorbed by the system. The first observation is that at sea states with significant wave height between 0.5 and 1 m, CL2 produces around 1 and 1.6 kW. For a generator nominal capacity of 30 kW this represents the region with more losses and thus the power production is little but still in the range of the numerical simulation results. The predictive controller results are presented in Fig. 31 linearised by the average of the mean electrical power productions of CL2 for each SS. The predictive controller presents an increase of electrical power in all the sea states between 13% and 65% in comparison with CL2. When taking into account the occurrence of each sea states, the weighted increase reaches 32%. The increase of power production can be explained by the fact that considerations are given to the entire PTO considering both the turbine and the generator during the optimisation. In Fig. 32 the efficiencies of the turbine, the generator and the total PTO are shown for each test in function of the pneumatic power. The turbine efficiency is higher in CL2 but focusing on the overall PTO efficiency, CL3 shows better performance. This validates the hypothesis raised in the analysis of numerical results that better consideration must be given to the generator behaviour. In addition, the efficiency of the PTO operating in CL2 during these low energetic wave conditions is around 30%, 35% for CL3. These scores are quite equivalent with the simulation results of Fig. 25, where the real sea states would correspond to those from SS1 to SS3. This means the numerical model is sufficiently well adjusted to predict the plant production under different environmental and control conditions.

Power quality
As for the numerical simulation case, an analysis of the quality of electrical power is done. Although only low energy sea states were experienced during the testing period, the values obtained allow a comparison of how the two controllers apprehend the power variations. Fig. 33 and Fig. 34 represent respectively the Pk2av power ratio and the standard deviation of the electrical power both with respect to the average of electrical power during all the tests.
Globally, it appears CL2 shows little benefits with lower values both in peaks and in standard deviation of the electrical power. Concerning CL3, results show power peaks a bit higher and slightly more variability in the power production than CL2. This study is quite limited due to the low values of power experienced, average electrical power of 25% of the nominal capacity at most. Thus it is difficult to draw conclusions on the impact of the predictive law on the power quality. As for the numerical results analysis, the predictive law captures more power in higher energetic group waves than the other advanced control law. There is clearly a lack of sea states and it should have been better to complete this analysis with higher energetic sea states to better understand how the predictive    law deals with peaks of power in a larger variety of sea states and compare them to the assumptions formulated in the numerical analysis.

Conclusion
The paper presented several scopes to analysis advanced speed control strategies applied to the PTO composed of the biradial turbine and the 30 kW induction generator installed in the Mutriku Wave power plant by providing both simulation results and sea trial ones. First the development of a W2W numerical model of the 9th chamber of the plant supported by operational data from the plant was more adapted to the realistic case and brought more accuracy to the simulation results.
Then it explained the development of 3 control laws for variable speed control, where 2 of them were state of the art and the latter was a predictive one. The main findings were that better PTO performance was associated to the predictive law over the two adaptive ones when focusing on electrical power production. When analysing in terms of the energy produced, CL2 produced 23% more than CL1 and CL3 produced 20% more than CL2. It appeared the good performance of the predictive law was due to the fact it was able to adapt to the incoming waves while being aware of system limitations (valve actuation, generator operation ranges) and that it optimised taking into account the overall PTO. Additionally, the lack of consideration of generator efficiency in the design of the controllers penalised the other laws. Also issues related to peaks of power were studied and the predictive law, although respecting the operation ranges, was more likely to operate at the components limits. One drawback of this kind of operation can be the accelerated wear of the electrical equipment and an oversizing of the power converter.
Real sea tests in the Mutriku Wave Power plant were performed for the best adaptive controller and the predictive one, being the first of its kind to be tested in a real environment to an OWC system. The tests lasted 2 weeks in the quiet wave conditions of June 2018. In the 7 sea states experienced during the testing period, the predictive law outran the adaptive one by 30%. The good performance of this CL offers promising results and could validate the implementation in a real control environment.
Future work will include a deeper analysis on the effect of including generator considerations in the adaptive control law and compare it with the predictive algorithm to quantify which is the dominant factor between the predictive feature and focus the control law on the generator performance. In addition, future development can be oriented to further investigate and tune the numerical model of Mutriku to better fit realistic conditions. It can include for example a comparison of the RAO obtained experimentally with the modelled one which will be insightful to deepen the understanding of the plant. Indeed, the MPC is as good as the numerical model. Also, the implementation into the control framework can be improved by implementing the controller directly into the PLC. Although the predictive controller could only be tested in low energetic sea states, the presented results are encouraging and justify to go on with a deeper analysis during a wider variety of sea states. The last testing phase within OPERA is to bring the best advanced control strategies tested in Mutriku and apply them to the Marmok A5 buoy installed in BiMEP, using the same PTO unit and in larger environmental conditions.