Prediction of Wave Energy Transformation Capability in Isolated Islands by Using the Monte Carlo Method

: In this work, a mathematical computer simulation model is used to predict the possible energy generated from different Waves Energy Converters (WECs) in the Canary Islands. The Monte Carlo Method is the computer simulation model proposed to predict the generated energy. The Waves Energy Converter systems analyzed in the study were, the Aqua Buoy, Wave Dragon and Pelamis converters. The models were implemented and validated, with the dataset of Gran Canaria deep water buoy. This buoy belongs to a network of buoys belonging to Spain’s State Ports and they cover a dataset period of 22 years. The research has concluded that it is possible to afﬁrm that the achieved model is a strong tool to compute the possible energy of any WECs, when the power matrix is known. The model based on the Monte Carlo simulation can be used in isolated islands of the Atlantic Ocean and can be extrapolated to other regions with the same characteristics.


Introduction
Spain has been a member of the European Union (EU) since January of 1986. The EU is an economic and political unification; currently there are 27 member states, including the UK, that are subject to the obligations and the privileges within the common judicial and legislative institutions. The EU recognize nine outermost regions, which are geographically far from the European continent. Only three countries (France, Portugal and Spain) have territory in the outermost regions, out of all the members of the EU [1,2]. These distant regions of Continental Europe represent a logistical challenge in many different fields, especially from an energy point of view. The energy production and transportation depend on imported fuels. In general, the electricity systems in the outermost regions are isolated. This fact increases the difficulty of its optimization [3][4][5][6].
The outermost regions in the North Atlantic Ocean are Azores and Madeira (Portugal). The Canary Islands (Spain) are the three Macaronesian archipelagos within the EU. All of them with volcanic origin and similar characteristics [7]. Figure 1 shows the Canary Islands that are formed by eight main volcanic islands located on the North Atlantic Ocean (27.5 • N-29.5 • N, 13 • W-18.5 • W), only 115 km away from the African continent. The population registered in the Canary archipelago at the end of 2019 was 2,153,389 inhabitants. Tourism is the economic driver of the Canary Islands.
Every single year millions of tourists visit the islands, except in 2020 due to the COVID-19 pandemic. The climate on the islands is generally warm and dry, creating arid almost desert-like conditions on the eastern islands (Lanzarote and Fuerteventura). Due to all these weather conditions, the Canarian archipelago has a serious water scarcity problem. Therefore, the service sector industry is only conceivable through the desalination of seawater, although this water is also used in agriculture and in supplying the population [3,4,[6][7][8][9][10][11][12]. Another difficulty in the archipelago is energy production. For the electrical network and electricity transportation, more than the 95% of the energy consumption comes from importations of fuels, with its respective problem of environmental contamination. The depths of the adjacent seawater are too great around most of the eight Canarian islands, which prevents the laying of submarine electrical cables between all the islands and with the African Continent, except between Lanzarote and Fuerteventura and Tenerife and La Gomera. The energy needs must be supplied locally in the majority of the islands [3][4][5][6]13]. For many decades the local government, in coordination with the Government of Spain, has sought an energy alternatives for remote islands. Due to this, it has focused a special interest on renewable energies (RE). Therefore, research and development institutions have been created to promote the use and generation of renewable energies. One example is the "Canary Islands Institute of Technology (ITC)" or "Gorona del Viento S.A" (www.goronadelviento.es, accessed on 30 August 2021), on El Hierro Island. This last institution has as a main objective to supply 100% of electrical energy from a hybrid renewable energy system (Wind-Pumped Hydro Power Station) in this small island [14].
Due to the geographical situation of the Canary Islands, it has abundant renewable sources such as wind and solar, which are increasingly used for the generation of electricity. In spite of several efforts that have been done to promote different kinds of the renewable energies (RE), little attention has been focused on wave energy, although different studies have shown that there is a considerable wave resource (25-30 kW/m) on the north and west coasts of the islands [15][16][17][18][19].
The wave energy will be one of the greatest ocean renewable energy resources, especially in islands and archipelagos where the bathymetric auspicious the generation of important wave energy concentrations [20][21][22][23].
Ocean energy is a fairly young sector and it is still developing. Over the coming years, one of these most powerful sources, wave energy, is estimated to play an important role in reducing CO 2 and greenhouse gas emissions worldwide, especially in the EU. Many EU coastal countries have an important ocean energy resource, which could help reduce the emission of polluting gases from the energy system and create new jobs. Currently in the Europe a wave energy resource has been identified in the range of 1000-1500 TWh [6,18,24,25].
The EREC (European Renewable Energy Council) considers that the marine energies, in particular wave energy, will be an important key in the goals of the European energetics strategy as it is relevant to the Directive 2009/28/EC on renewable energy. The EU has committed to reducing its polluting gas emissions by around 80-95% by 2050 [6,[25][26][27].
The wave energy sector demands attention respect to other coastal energy technologies, for example offshore wind energy, and there are many steps in this direction [28]. Stratigaki [29] discuss one of the most important project about wave energy in Europe, the WECANet (Wave Energy Converter Array Network) European COST Action (www.wecanet.eu, accessed on 30 August 2021). The project is supported by the European COST Association (COST stands for "European cooperation in science and technology") and was introduced in 2018. WECANet is a consortium of professionals from 30 European member countries and United Stated on America. The aim of this project is to create a solid network platform focused on current challenges, identifying barriers and helping the development of wave energy research and industry.
Numerous studies and assessments about the wave energy potential have been implemented in the Macaronsian region, included the Canary Islands. On the one hand, the computer third-generation (3G) wind-wave models, WAve Modelling (WAM) and Simulating WAves Nearshore (SWAN) [15,17,[30][31][32][33], has been used widely as an essential tool for the development of these studies. On the other hand, other numerical systems to wave forecasting in the western Atlantic Ocean based on soft computing techniques have also been used [6,34].
As general knowledge into the renewable energy industry, wave power converters (WECs) designs have not been developed for being used on a commercial scale. Moreover, an extensive number of WECs demonstration prototypes have been carried out in the EU, almost large-scale prototypes, but they ran into technical or financial problems. Despite the difficulties, the wave sector is not stopped by the setbacks and continues growing with the economic help of European governments up to the time when a competitive commercial prototype is achieved [29,35,36].
In addition, an important part in the progress and expansion of RE based on waves is the prediction of the behavior of any Wave Power Converters (WECs) based in their power matrix, especially in island regions with high wave potentials, external fuel dependency and isolated electricity systems. For that motive, the main purpose of this study is to determine the behavior of the Monte Carlo simulation to predict the energy generated from different Waves Energy Converters (WECs) in the Canary Islands, in order to establish the advantages entailed by the use of this mathematical method, in the prediction of the possible wave energy generated from a WEC. The models were carry out with the dataset from Gran Canaria buoy (2442). This buoy belongs to network buoys of the Spain's State Ports and covers a dataset period from 1997 to 2019.
Monte Carlo method is a simulation that relies on recurrent random sampling and statistical analysis to obtain numerical results. The essential concept is to use randomness in the experiments, for which the particular result is not acknowledged in advance [37]. The mathematical simulation model can be used in many fields of engineering such as, energy model simulation [38,39], renewable energy generation, prediction and economic analysis (Wind, Photovoltaic, Wave Energy, etc.) [40][41][42][43][44][45][46][47][48][49]. The computational simulation is really useful in ocean industry too, being used in different studies such as, marine and ships structures [50][51][52][53], breakwater analysis [54], offshore wind turbine [55], maintenance operations of Wave Energy Farm [56], etc.
Other applications in the maritime field can be seen in the proposition of a new ships route decision-making strategy [57,58], reducing the accident risk in marine transportation [59], ships engine emission control [60], offshore transportation [61], ship fire consequence, etc. [62].
Some research area conditions such as regulatory frameworks, environmental impacts, social dimension, fishing activity areas, protected areas, maritime routes, military zones, bathymetry, submarine cables, pipeline protection areas, have not been considered in this study.
The goal of this paper is intended to predicting the energy transformation capability, based on historic wave data and considering the random nature of this behavior through a combination of bivariate Weibull distributions and Monte Carlo Method.
This work is structured in eight sections. After the introduction, the wave data sources from the buoy and the WECs selected for the study are shown in the section two and three. The fourth and fifth sections display the general modeling procedure of the Weibull distribution. Monte Carlo simulation of the power conversion are shown in section six. Finally, the validation the proposed model and conclusions of the whole research are given in the sections seven and eight, respectively.

Data Sources
Currently in the Canary Islands there are only five buoys that provide data on the sea waves in real time, two deep waters buoys and three nearshore buoys. The deep waters buoys are: Gran Canaria buoy (2442)  Data from the Gran Canaria buoy (2442) was taken to predict the energy generated from different WECs with Monte Carlo simulation in this study, because it was the buoy with the best position compared to the other buoys. It is possible to carry out this statement because the north and west of the Canary Islands are the most energetic maritime zone, due to the influence of the trade winds. Buoy 2442 has another advantage, it is the least affected by the shadow of the islands. The buoy is placed in the Northwest of Canary Islands, located at 28.20 • N and 15.78 • W. The geographic location of the buoy is displayed in Figure 1 [6,15,63].
The Gran Canaria buoy belongs to the deep water buoys network of the State ports of the government in Spain. This buoy is mooring at more than 200 m of deep waters. Note that, the buoy 2442 is not affected by refraction or shoaling due to the lack of continental shelf of the Canary Archipelago [6,15,63].
The data set used covers a historical period of 22 years, from 1997 to 2019, when the Gran Canaria buoy was removed to harbor for maintenance. The whole dataset was divided into a training set, which comprises data corresponding from 1997 up to year 2016, and a validation set, including data from year 2017 to year 2019.
The measurements of variables as, mean peak period, T, and wave significant spectral height, H, were carried out every three hours in the Gran Canaria buoy. The measurement accuracy was ±0.05 s for T and ±0.05 m for H [63]. This study has only been possible thanks to the collaboration received from the Harbors of the State of Spain, which provided the entire data set.

WECs Analyzed
The combination of wave energy systems with other renewable energies systems is particularly interesting for isolated electricity systems, which in many instances are totally isolated from continental electricity systems [64].
A retrospective of the wave energy industry shows an extensive variety of WECs concepts in development in dissimilar countries around the world [64,65]. Ahamed, McKee and Howard [65] showed that until 2019 there were more than one hundred wave energy experimental projects. All these projects have been installed and tested in dissimilar countries around the world, for example, USA, Australia, China and New Zealand. In Europe the counties involved in waves projects are Sweden, Portugal, Italy, Norway, France, Spain and UK [65,66].
Selecting the optimal deployment places involves previously feasibility studies and mathematical simulations to support the decision [67,68]. In this work, the Monte Carlo Method has been used to determine the behavior in deep water near to Canary of different tested WECs prototypes.
The number of devices considered on work are three: Pelamis (750 kW), Wave Dragon (7000 kW) and Aqua Buoy (250 kW), with different dimensions and power take off (PTO). The operating principle used for the wave energy transformation systems to change the sea wave oscillating motion into electrical energy are two, the hydraulic motor in the Pelamis case and hydro turbine, in the cases of Wave Dragon and Aqua Buoy [69][70][71][72].
The power output of any WECs in different sea states, taking into account the mean peak period, T, and the significant spectral height of the waves, H, are called power matrices. Figure

Bivariate Weibull Distribution Fitting
For describing the stochastic behavior of the wave, given by the numerical values of the significant spectral height, H, and the mean peak period, T, a bivariate Weibull probability distribution function is fitted: where c 1 , k 1 , c 2 , k 2 , and c 12 are the distribution parameters, and B 0 is the zero-order Bessel function. The corresponding parameters are obtained, through the Simultaneous Maximum Likelihood Estimate method [30,63,73,74], for each week in the year. Points in Figure 3 show the obtained values for each week, of the training set, from 1997 up to year 2016. As can be seen, all the coefficient values move around certain average. Coefficients c 1 and c 2 , which represents the values of significant spectral height and mean peak period, respectively, where the maximum value of probability distribution is reached, show a noticeable season variability, with higher values from October to March. Seasonal variations of coefficients k 1 and k 2 , which determine the shape of the probability distributions, show a light increment from May to September. Finally, it should be also noted that coefficient c 12 does not evidence seasonal variations, indicting no changes in the interaction between significant spectral height and mean peak period.

Bivariate Weibull Coefficient Models
In order to calculate the parameters of the bivariate Weibull coefficients for each day computed from the year start, d, data from 1997 to 2016 were used. The relationship was expressed by a sixth-order polynomial equations, in the form: Coefficients for these models, β 6 , . . . , β 2 , β 1 , β 0 , were obtained though stepwise least square method, with backward elimination where the model starts with all the involved variables followed by an iterative deletion of the less statistically significant variable. The process is repeated until no further variable can be deleted without a statistically nonsignificant loss of fit. Equations used for fitting the models and obtaining the corresponding standard deviation are shown in Appendix A.
Obtained models were summarized in Table 1. As can be seen, all the probability associated with the Fisher test, p(F), in the analysis of variance (ANOVA), are lower than 0.01, meaning that there is a statistically significant relationship between the variables, at the 99% confidence level. The maximum values of the probability associated to the corresponding t-Student tests, show that all the coefficients kept in the models were statistically significant at 99% confidence level. Finally, the values of coefficient of determination (R 2 ) and mean absolute error (MAE), indicate that the models as fitted explain only a low percent of the variability in the coefficients. It must be pointed out that a low R 2 value does not automatically imply a bad model, especially for fields of study having an inherently greater amount of unexplainable variation (such climate phenomena [75,76]). In these cases, models can be accepted even with a low R 2 , given that independent variables are statistically significant. However, the high variability in the obtained models is the main reason for using the Monte Carlo method, which allows us to evaluate the random behavior of the modeled variables.  Figure 3 shows the predicted values (thick lines) and the 95% confidence intervals (thin lines). It can be noted that the predicted values of c 1 (which determines the value of significant spectral height were the probability function reaches its maximum) moves from 1.5 to 2.1, with highest values in winter (November-February) and lowest values at the end of the summer (August-September).
The predicted values of k 1 range from 4.1 to 6.3, which determines a bell-shaped probability function, on the significant spectral height dimension. This model reaches its maximum in June-July and minima in February and November.
The values of c 2 determines the value of the mean peak period where the probability function achieves its maximum. The corresponding predicted values of c 2 range from 8.1 to 12.6. The highest values for this coefficient take place in winter (November to February), while its lowest values, in summer (June and July).
The predicted values of k 2 range from 5.5 to 8.8. Two minima can be observed, for these predicted values, corresponding to March-April and October-November, and two maxima, in December-January and June-July.
The last model (c 12 ) is a special case, because all the terms, except the constant, are statistically non-significant. Consequently, the equation is just a constant, meaning that there is not a significant variation of c 12 with the time and all the observed variability has a random nature.
All the fitted models, and the corresponding standard deviations of their residuals, are shown in Appendix B.

Random Values Generation
In order to perform the Monte Carlo based simulations, random generators should be defined for each bivariate Weibull distribution coefficient. These generators consist on normally distributed random numbers, with means and standard deviations as computed, in the previous subsection, for each coefficient. This can be formalized through the equation: whereĉ is a random value of coefficient c, N (•, •) is the normally distributed random number, and µ(c) and σ(c) are the estimated mean and standard deviation of coefficient c. Figure 4 shows thirty random generated values, for each bivariate Weibull, in the whole interval of days of the year. As can be noted, there is a good matching between the observed (see Figure 3) and predicted values for all the coefficients.

Power Conversion Computing
The power conversion capabilities of the analyzed point, was studied with three different devices: Aqua Buoy, Wave Dragon and Pelamis. Each of them has the corresponding power matrix, which are shown in Figure 2.
With every randomly generated set of bivariate Weibull coefficients, the cumulative probability, in each interval of the matrix (H l ≤ H ≤ H u , T l ≤ T ≤ T u ), is determined, by the equation: where f (H, T) is the Weibull probability distribution function, given by Equation (1). This probability is numerically computed, by using the fifth order Gauss-Legendre quadrature for double integrals: where H i and T j can be computed by the expressions: and w i and ξ i are the corresponding weights and evaluation points (see Table 2) [77]. The total converted power, P, for each device is then computed through the expression: where m and n are the intervals for significant spectral height and mean peak period; φ represents the conversion power for the considered interval and [H i l , H i u ] and [T i l , T i u ] are these intervals. Figure 5 shows the graphical representation of the computed converted power, for the 30 random bivariate Weibull distributions, through the year, and the corresponding mean value and 90% prediction interval. As it can be appreciated in Figure 5, the conversion of sea wave oscillating motion into electrical energy, is not excessively high in the selected point for the study, if the nominal power of the Pelamis (750 kW), Wave Dragon (7000 kW) and Aqua Buoy (250 kW) WECs are taken into consideration.

Validation
In order to validate the proposed model, the prediction interval for each month was computed for each month of the year. The obtained outcomes were compared with the actual computed values for years from 2017 to 2019 (see Figure 6). It is worth mentioning that the entire observed values fell into the predicted intervals, except those corresponding to January and February 2018 in the Pelamis and Aqua Buoy converters, respectively. However, these values are very close to the upper limits.
An interesting issue arises from comparing the performance (i.e., the ration between the actual and nominal converted energy) of the three evaluated converters. In spite of the remarkable contrasts in the amount of converted energy, there are not noticeable differences in their performances. As can be seen (Figure 7), te average converted power ranges from 6% to 18%, while the maximum converted power ranges from 12% to 30%. The cause of this similarity if the divergence between the maximum probability values for the waves significant spectral wave and mean peak period (from 1 to 3 m and from 6 to 15 s, respectively), and the values were the maximum energy conversion take place for the three converters (see Figure 2).

Conclusions
The transformation of sea waves oscillating motion into electrical energy, signifies a complex and expensive challenge. Currently there are abundant WECs concepts, but there is no consensus about a unique design approach. It might be caused by the high cost of the obtained kilowatts.
This research presents a mathematical simulation based in Monte Carlo Method, to predict the energy generated from different WECs in Canary Islands. The model was implemented and validated with dataset from Gran Canaria deep waters buoy (2442), which belongs to the network buoys of Spain's State Ports.
In the bivariate Weibull distribution fitting used in the study, the distribution parameters c 1 , k 1 , c 2 , k 2 , and c 12 , moves around certain average, a season variability can also be detected, in c 1 and c 2 parameters.
In the polynomial regressions used to model the relationship between the coefficients values and day number, it is possible to see that the in the regression model for the coefficients c 1 , k 1 , c 2 and k 2 , the probability value of the ANOVA table was lower than 0.01. This means that there is a statistically significant relationship between the variables at the 99% confidence level. With respect to c 12 , it is possible to affirm that there not a significant variation with the time of the coefficient and all the observed variability has a random nature. It should be recalled that the small R 2 values obtained in the study does not imply a wrong model, particularly for very noisy models, such as those based on climatic observations. Before calculating the wave power conversion of each WEC, it is necessary to perform the Monte Carlo-based simulations for each bivariate Weibull distribution coefficient, which simulation shows a good matching between the observed and predicted values for all the coefficients.
The power conversion capabilities or power take off (PTO) of the study devices, Pelamis, Wave Dragon and Aqua Buoy on the analyzed point (buoy 2442) was calculated with the power matrices of each one. The transformation of ocean wave oscillating motion into electrical energy, is not so high in the selected point for the study, taking into consideration the nominal power of the WECs. This results do not differ from the quality of the model based in the Monte Carlo Method, being demonstrated its effectiveness in the validation process.
The proposed approach has two main advantages. In the first place, considering the variability in the wave behavior allows to predict not only the expected converted energy but also its actual confidence interval. This plays a key role in the design of power systems were minimum delivered energy should be taken into account. In the second place, this methodology offers a model for the wave behavior, based on the bivariate Weibull distribution, and uses this model for simulating the performance of a given wave energy converter. Separating both steps, it is possible to evaluate new converters without fitting a new model, which cannot be done through other approaches, such as time series-based modeling of the converted power, directly evaluated from the raw measured wave data.
The model developed in this study could be an effective tool to calculate effectively and precisely the PTO of different WECs at any position in deep ocean waters, preferably on isolated islands of the Atlantic Ocean. The model will be useful in forecasting renewable energies, this energy can be introduced into the electrical network, used directly in desalination plants or in combination with other kind of renewable energies.

Acknowledgments:
The authors would like to thank Puertos del Estado (Spain's State Ports) for having provided the datasets of the buoys used in the study and the Technological Institute of the Canary Island (ITC) for the given support.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. General Modeling Procedure
For modeling the relationship between the coefficients values and day number, measured from the year start, d, polynomial regressions, in the form: were fitted. Coefficients β = {β 6 , . . . , β 2 , β 1 , β 0 } T were obtained by using the expression: where the days (independent variable) values, of the training dataset, are arranged into the matrix, being: and n and p are, respectively, the number of training dataset samples and the number of considered model coefficients.