A Dynamic Stochastic Hybrid Model to Represent Significant Wave Height and Wave Period for Marine Energy Representation

This paper presents a methodology to represent ocean wave power generation based on real data observation for significant wave height (SWH or Hs) and wave period (WP or T). This technique is based on a hybrid model, which considers Fourier series and stochastic differential equations, allowing a continuous time representation of the random changes in the parameters associated with wave power generation (Hs and T). The methodology is explained, including estimation methods and a validation procedure. The data series generated by the models erre used to create simulated wave power output applying a transformed matrix and a theoretical model. The results validate the utilization of this technique, when the objective is to obtain a robust dynamic representation of a random process, oriented to linear studies.


Introduction
The inclusion of wave energy in electrical power systems (EPS) is a current research topic at global level, due to its energy potential-the prospective world production has been theoretically estimated at more than 100,000 TWh per year [1]-and the increasing interest to include renewable energy supplies [2]. For example, North American, European and Asian countries have installed wave energy conversion facilities [3], e.g. Ireland with the ocean energy buoy [4]; the United Kingdom with Pelamis, Ceto I and Ceto II, Oyster and Limpet [5,6]; the United States with PowerBuoy [7]; and Denmark with Wave dragon and wave star [6,8]. Additionally, in African and Oceanian regions, some studies related to implementation and integration of wave energy conversion are being developed [9]. Even though there is increased interest in this energy source, there still exists geographic zones with significant potential, which have not been intervened, for example, the case of Chile in South America [10], which has access to the sea in its whole continental territory.
In all those implemented examples and in future projects, it is necessary to study before a real implementation, the effect of integrating wave energy into the grid, because this type of power induces random perturbations that can trigger damage effects in the system at operational level. Due to the

Problem Formulation
This paper proposes robust dynamic stochastic models of H s and T, which allow representing wave power. These models are oriented to linear studies of inclusion, such as small signal stability analysis, where the linearized equation of the complete system are investigated with the objective to maintain the stability at a generation level. Thus, the models to be presented consider simplified components, i.e., they consider linear and stationary elements, despite the non-stationarity and non-linearity of the ocean behavior.

Theoretical Model
The proposed models characterize the variation of H s and T over time considering real data observation for the parametric adjustment and validation. The methodology consider two parts: a first analysis based on Fourier series, which is used to describe the seasonal behavior of waves in deep waters [23], where the theory of linear waves can be applied [13]; and a second part corresponding to the representation of the randomness inherent to the oceans, which is characterized by the Ornstein-Uhlenbeck (O-U) process. Mathematically, this is expressed in Equation (1) for H s and in Equation (2) for T.
where H t is the H s at the instant t, T t is the T at the instant t, f i is frequency associated to seasonal changes, U t O-U process represents the noise, and α, γ i , δ i , β, ζ i and κ i are parameters to be estimated.

Numerical Model
To compute numerical values, according to our theoretical model, Equations (1) and (2) can be written in their differential form using Ito's formula [24], and are expressed in Equations (3) and (4), respectively.
The first stage of the proposed procedure consists of decomposing the data sample into its frequency spectrum f i , using the periodogram, which corresponds to an approximation of the spectral density of the data series. With those frequencies, it is achievable to establish an equation composed of a sum of sine and cosine, linear in the parameters α, γ i and δ i in Equation (1), and in the coefficients β, ζ i and κ i in Equation (2), allowing the use of linear regression techniques for parametric estimation as least square. Only the significant frequencies must be considered, which correspond to those frequencies with a non-zero associated parameter in the sinusoidal function.
The second stage consists in associating SDE to the residuals of the Fourier adjustment, which os the Ornstein-Uhlenbeck process. In general, SDE can be used to simulate the continuous behavior of a variable, which is why in engineering topics it is common to find applications to power systems for short-term studies and real time studies [25][26][27], such as the case of small signal stability [28]. However, as an application to the wave generation area, it is not a highly exploited resource [29].
The SDE can be considered as a generalized case of an ordinary differential equation, where the variable under study, state or dependent presents a random behavior in its evolution. One of the most common SDEs used to describe physical phenomena is the so-called Ornstein-Uhlenbeck (O-U) process, with Equation (5).
where U t corresponds to the O-U process, which represents the state variable or dependent, η < 0 and ν > 0 are parameters and W t is the Brownian motion or Wiener process. Besides, given the initial condition for the process U t , it follows a normal distribution, i.e., ).

Model Validation
To adjust the parameter of the O-U process and validate it, the residuals left by the Fourier analysis are used as data. To validate the results once the model is adjusted, the method called uniform residuals analysis [30] is used, which in simple words determines if the actual observations of a certain phenomenon are the result of a process with the characteristics of the proposed model.
With the parameters of the models of H s and T estimated and validated, the SDE of the models are simulated. The data generated allow a representation of power generation, transforming the data with physical and empirical models.

Summary
To summarize, the methodology is described in the following stages: • Have a data source of H s and T.

•
Apply the Fourier analysis to the data sample, i.e., estimate the significant frequencies by means of the periodogram and the amplitudes associated to those frequencies, using the least square estimation technique.

•
Adjust the residuals of the previous estimation to an O-U process, using them as observation for the parametric adjustment, which is made by quasi-maximum likelihood technique [31] for the drift parameter η, and quadratic variation for the diffusion coefficient ν. • Validate the resulting model with the techniques of uniform residuals [30].

•
Simulate data series of H s and T, using the Milstein scheme for the corresponding SDE [24]. • Transform the simulated series into power generation information by means of the physic model for generated power in deep water [22] and the power matrix of the floating oscillating water column device [16].

Implementation to the Chilean Case
This section presents the systematic application of the modeling methodology. The implementation allowed the generation of future data series of H s and T, which provide information on the generation power through the transformation of these series of data.

Data Sources and Characteristics
The series of H s and T correspond to H s measurements in meters and T in seconds, with a resolution of 3 h, which were generated by two buoys located on the coast of Valparaiso (see Figure 1, red marks) belonging to the Servicio Hidrográfico y Oceanográfico de la Armada de Chile (SHOA). The dataset was used with the objective of illustrating the application of the methodology. There are no better samples of data for the Chilean coast. Table 1 summarizes the specifications of this information:

Parameter Adjustment
According to the methodology presented in the Section 2, the frequencies of Equations (1) and (2) must be estimated by means of the periodogram, for each sample of data. The significant frequencies (Hz) estimated are shown in Tables A1 and A2 in the Appendix A. Then, with that information, it was possible to estimate the coefficients α, γ i and δ i in Equation (1), and the coefficients β, ζ i and κ i in Equation (2) for each sample of data, using the least Square technique. The estimated values for each parameter of each model are summarized in Tables A3 and A4 in the Appendix A.
The residuals left by this adjustment were fitted by the O-U process in Equation (5), hybridizing the methods of Fourier analysis and SDE. These residues corresponded to the observations of the continuous stochastic process O-U, which was used to estimate the parameters of the proposed process as a model.
For the parametric estimation of η, the quasi-maximum likelihood technique [31] was used. Under this method, theη estimator is given by Equation (6).
where U i∆ , i = 1, ..., n are equidistant observations of the O-U process and ∆ is the time interval between observations. It was not possible to estimate parameter ν using this method, and in this case the estimation was done by means of the quadratic variation of the process. The estimator was obtained by equating the theoretical expression of the quadratic variation with its discretization. For the O-U process, the estimatorν is given by Equation (7).
where t corresponds to the total observation time interval.
The values of the parameters delivered by Equations (6) and (7) for each model and each respective sample are presented in Tables 2 and 3.

Validation Procedure
Once the parameters were estimated, the model Was validated. The method used for the present work corresponded to the analysis of uniform residuals [30]. This technique determined whether the data being adjusted could actually be described by the proposed process. The result of this analysis was graphically observed by means of a quantile-quantile plot, which compared the experimental distribution of the data with the theoretical distribution of the stochastic process. If the identity line were obtained, then it could be concluded that the model satisfactorily adjusts the observed data. Figure 2a-d shows the quantile-quantile plot of this analysis for each sample data and each model.
The quantile-quantile plots in Figure 2a,c correspond to the settings for Sample 1, which consisted of 248 observations. It was observed that the points in the graphics followed approximately the identity line, thus it was affirmed that there is no reason to doubt these adjustments. For the plot in Figure 2b,d, corresponding to Sample 2, which consisted of a total of 3207 observations, the line was not precisely the identity line. The deviation of the curve suggests that the proposed process failed to describe certain characteristics of the data, removing the precision of the model. The difference between the results of Samples 1 and 2 might be due to the fact that they were very different in size: Sample 2 consisted of many data. Tthe nonlinear and non-stationary behaviors of the oceans might be very strong and do not allow the correct fit of the proposed hybrid model. This was not observed in Sample 1, which corresponded to a much shorter observation interval in which it is possible that nonlinear and non-stationary effects were not fully captured.

Results
Once the model was validated, the projections of present and future values of H s and T must be done to apply the respective transformations to the series and obtain simulations of generated power.

Simulated Series of H s and T
To generate the trajectories made by Equations (1) and (2), the numerical scheme of Milstein was used, corresponding to an approximation up to the second-order of the Ito-Taylor expansion of the SDEs under study in Equations (3) and (4) [24].
Equations (8) and (9) express the simulation schemes for the H s and T models, respectively.
where W t j+1 − W t j corresponds to the increments of the Brownian motion, which is simulated by means of the relation W t j+1 − W t j = (t j+1 − t j )N j , where N j ∼ N(0, 1). Both equations were simulated with a time interval (t j+1 − t j ) = 0.12 h, allowing the generation of more refined series of H s and T than the originals. In total, 500 simulations were performed for both equations, as plotted in Figure 3a-d. By way of comparison, the average of the 500 simulated trajectories is plotted in black together with the graph of the real data trajectory in red. The blue line in the graphs indicates to the left side the 80% of data used for the parametric adjustments and to the right side the 20% of data used for comparing the performance of the model at the moment of forecasting future values. This allowed measuring the prediction errors of the settings in the next stage.  In the case of T, Figure 3c, corresponding to the shorter sample, shows that the variation-understood as amplitude around an imaginary moving average-of the model was much lower than that of the real data; this effect can be seen more clearly in Figure 3d. However, the ability of the model to follow the trend of the wave period stands out.

Hybrid Model Performance
The models were subject to a validation process and the errors of the results given by the models were measured by the mean percentage error (MPE). Table 4 summarizes the values obtained for the MPE of the H s and T models for both data samples.  Table 4 shows the adjustment error, corresponding to the comparison of the 500 trajectories generated with the 80% of total data used for the parametric adjustment. On the other hand, the forecast error corresponds to the comparison with the 20% of total data used for the respective part of the average of the 500 trajectories.
The adjustment errors of H s and T models of both data samples were acceptably low: less than 3% for the adjustment and less than 9% for the future forecast. However, these models use linear wave theory, and it is expected that the error of forecasting increases due to the non-linearities associated to the physical measurements, but, as the paper seeks a robust model that can be useful in small signal analysis, the measurements errors are in an acceptable range.

Physical Power
The interest variable of power generation was obtained by means of two methods. The first corresponded to a transformation of the data generated in the previous section through the physical relationship for the deep-water wave-power level [22], expressed in Equation (10), where it takes as an input of H s (SWH) and T (WP) the series generated by the proposed models.
where ρ = 1025 kg/m 3 is the mass density of sea water and g = 9.8 m/s 2 gravity acceleration. The series for the power generated under Equation (10) for both samples of data are presented graphically in Figure 4a,b, respectively.
As shown in Figure 4a, a good agreement of the present data was observed, following the slope changes correctly. For future projection, although the change in trend was well described, the random variation was greater in the real data than in the simulated average trajectory.
On the other hand, Figure 4b shows a remarkable decrease in the description of variability by the model, in both present and future values. Regarding the characterization of trend change, a good performance was observed.  Graphically, a decrease of the variance should be noted. This was due to the fact that the simulated observed trajectory corresponded to an average of 500 simulations and not to a single trajectory.

Floating Oscillating Water Column
The second method to obtain a series representing the power generation corresponded to a transformation of the data generated for H s and T through the power matrix of the floating oscillating water column device (OWC) [16], which was obtained experimentally together with the Pierson-Moskowitz spectrum analysis. The selection of this device was based on the fact that the operating ranges of the machine matched the data provided by the buoys. The power matrix for this device is presented in Table 5. In Figure 5a,b, the black path corresponds to the power trajectory generated by the average of the 500 simulated H s trajectories. In red, the power generated by the actual data is plotted and the blue vertical line represents the set point to the left, and the forecast to the right.
The graphical results are analogous to those observed in Figure 4a,b.

Conclusions
The traditional techniques that are used to capture randomness in EPS consider deterministic models mainly. The results or estimations are validated with Monte Carlo-like simulations. In our case, we use stochastic differential equations to make sure we capture the randomness that is self-sustained in time and inherent to marine resources. Here, the validation process for the obtained parameters is based on the use of the residual method for these particular differential equations. This should be regarded as very important because it allows having a robust method to assess quality in the estimated parameters.
The objective of this work is to propose a stochastic dynamical model, which is based on Fourier analysis, applicable to linear wave theory, and stochastic differential equations. With this approach, it is intended to represent the seasonal changes in the waves and their stochastic characteristics. The proposed models were applied to the characterization of two fundamental variables in the wave energy generation: significant wave height (SWH or H s ) and wave period (WP or T). Since we had real observations of these two variables, they were used to adjust the parameters involved in the models by means of statistical techniques to obtain a characterization of the particular region in which the models were applied.
Wave power was obtained using the data series generated by the proposed models in two different procedures: a theoretical equation valid for deep water and the power matrix of the floating oscillating wave converter (OWC), which was chosen given its conditions of operation, applicable to the region where the real measures of H s and T were obtained. The order followed in the procedure to generate series of power data facilitated the use of linear models, since it avoided directly modeling the power, which has a non-linear behavior.
The paper includes the systematic explanation of the methodology described above, explaining the methods of parametric estimation, validation of models, simulation of present and future values, and adjustment of error measures, which is done by means of the mean percent error.
The results show the good performance of the proposed hybrid models for H s and T when adjusting present values in time windows in the range of one month, with an error of 1%. For larger samples, the results show a failure in the precision of the method since it could not capture nonlinear and non-stationary components of the oceans, presenting an error close to 3.6%. The same is true for projections of future values, increasing the error to 8%. It can be concluded that a linear approximation of the ocean behavior is useful and valid to be applied to inclusion studies of this energy source in the short term, such as in the case of small signal analysis for power systems, where the equations are reduced to the linear case. Within this same context, an additional advantage of the proposed methodology, based on techniques in continuous time, is to generate data series with temporal resolutions finer than the series of data delivered by the measurement devices.