Wave Excitation Force Estimation of Wave Energy Floats Using Extended Kalman Filters

TL;DR

This paper explores the use of Extended Kalman Filters to estimate wave excitation forces on wave energy converters from measurement data, enhancing control strategies for wave energy devices.

Contribution

It introduces and compares two novel EKF-based methods for estimating excitation forces, including a direct state inclusion and a harmonic oscillator model, validated with experimental data.

Findings

Both methods successfully estimate excitation forces under various wave conditions.

The harmonic oscillator model shows robustness to measurement noise.

Sensitivity analyses reveal the impact of model parameters and sampling rates.

Abstract

In many advanced control strategies the wave excitation force is key to determining the control input. However, it is often difficult to measure the excitation force on a Wave Energy Converter (WEC). The use of Kalman filters to estimate the wave excitation force based on readily available measurement data can potentially fill the gap between the development of WEC control strategies and the data that is available. Two different estimation methods using an nonlinear Extended Kalman Filter are tested on experimental wave tank data for a heaving semi-submerged float. The first method relies on directly including the excitation force as a state in the first order dynamics---which allows the "random walk" of the Kalman filter to identify an estimate of the excitation force. The second method of estimation involves modeling the wave excitation force as a harmonic oscillator comprised of sinusoidal components. Both methods are evaluated for a variety of incident waves and additional sensitivity analyses are performed to investigate the susceptibility of these estimation methods to changes in the model, measurement noise, and sampling rate.

Figures

Tables

Table 1. Table 1: List of experimental runs from [ 15 ] that are considered in this work with corresponding significant wave height, H s subscript 𝐻 𝑠 H_{s} , and peak period, T p subscript 𝑇 𝑝 T_{p}

Run No.	$H_s$ mm	$T_p$ [s]
1	175	1.42
2	50	1.90
3	175	1.90
4	300	1.90
5	50	2.37
6	175	2.37
7	300	2.37
8	50	2.69
9	175	2.69
10	300	2.69
11	50	3.16
12	175	3.16
13	175	3.95
14	175	4.74
15	175	5.53
16	375	2.06
17	375	2.37

Equations

$$M_v\ddot{y}$$

$$F_M=-K_{mu}y,$$

$$F_B=-\rho gAy,$$

$$F_D=-\frac{1}{2}\rho A_p{C}_{D}V|V|,$$

$${\dot{X}}_r(t)=$$

$$F_R=$$

$$F_{ex}(t)={\int }_{-\infty }^{\infty }{F}_{IRF}(\tau )y_w(t-\tau )d\tau .$$

$$\dot{\text{x}}(t)$$

$$\tilde{\text{y}}(t)$$

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ x_r\end{array}\right]=\left[\begin{array}{c}{x}_2\\ \frac{1}{M_v}F(x_1,x_2,F_{ex},\dot{y},C_r{x}_r)\\ A_r{x}_r+B_r\dot{y}\end{array}\right]$$

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\\ {\dot{x}}_r\end{array}\right]=\left[\begin{array}{c}{x}_2\\ \frac{1}{M_v}F(x_1,x_2,x_3,\dot{y},C_r{x}_r)\\ 0\\ A_r{x}_r+B_r\dot{y}\end{array}\right]$$

$$F_{ex}=\mathbb{R}\left[\frac{H}{2}{F}_X({\omega }_r)e^{i({\omega }_{r}t)}\right],$$

$$F_{ex}=\mathbb{R}\left[{\int }_0^{\infty }\sqrt{2\Delta \omega S(\omega )}{F}_X(\omega )e^{i(\omega t+\varphi (\omega ))}d\omega \right],$$

$$f_{ex}$$

$${\ddot{f}}_{ex}$$

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{f}}_{ex}\\ {\ddot{f}}_{ex}\\ {\dot{x}}_r\end{array}\right]=\left[\begin{array}{c}{x}_2\\ \frac{1}{M_v}F(x_1,x_2,f_{ex},\dot{y},C_r{x}_r)\\ {\dot{f}}_{ex}\\ -{\omega }^2{\dot{f}}_{ex}\\ A_r{x}_r+B_r\dot{y}\end{array}\right],$$

$$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{f}}_{ex}\\ {\ddot{f}}_{ex}\\ {\dot{x}}_r\end{array}\right]=\left[\begin{array}{c}{x}_2\\ \frac{1}{M_v}F(x_1,x_2,\Sigma {\bar{f}}_{ex},\dot{y},C_r{x}_r)\\ I{\dot{\bar{f}}}_{ex}\\ -{\Omega }^2{\dot{\bar{f}}}_{ex}\\ A_r{x}_r+B_r\dot{y}\end{array}\right]$$

$$y(t)=\frac{F_{ex}(t)}{||F_X({\omega }_p)|{|}_2}$$

$$NMSE=1-\frac{||s(t)-\tilde{s}(t)|{|}_2}{||s(t)-\bar{s}(t)|{|}_2}$$

Full text

††journal: Journal of Ocean Engineering

1 Introduction

Wave Energy Converters (WECs) are devices that generate electricity from ocean waves. Some non-exaustive reviews of the widely varying devices and modeling methods for WECs can be found in [1, 2, 3, 4].

This work will focus specifically on point absorbing wave energy converters. Point absorbing type WECs rely on a small body, when compared to the incident wavelength, to be actuated primarily in a single degree of freedom to achieve power generation. A point absorbing WEC was chosen for this work as point absorbers typically offer significant deployment versatility and simplified heave dynamics.

Modern control strategies rely on state information to calculate an optimal control input. As addressed in [5] and [6], this is possible through the use of Kalman filters as optimal state observers. Kalman filters can also be used to estimate the wave excitation force which is necessary for many controllers [6, 7]. Even though Kalman filters are commonly used to estimate the states for control systems, the implementation of a Kalman filter on a WEC is not trivial. This work addresses some of the current barriers to state observation and excitation force estimation of wave energy converters using an Extended Kalman Filter (EKF), which is the nonlinear adaptation of the Kalman filter.

The use of an EKF to determine wave excitation force using external pressure sensors is shown in [8]. However, the implementation of additional external pressure sensors may be impractical for many commercial WEC applications. Alternatively, [6] shows a way to use a Kalman filter to estimate the exitation force by modeling the excitation force as a disturbance in the WEC equations of motion. The disturbance used in [6] describes the excitation force as a sum of sinusoids as is typical in WEC modeling [4]. This method of estimating the wave excitation force will be referred to as the disturbance estimation method. Alternatively, a linear Kalman filter is implemented in [7] where the excitation force is identified by including the force as a state with no dynamics and thereby allowing the “random walk” of the Kalman filter to identify the best estimate of the state directly. This methodology of implementing the Kalman filter to identify the wave excitation force will be referred to as the direct estimation method. The work of [7] also develops a receding horizon estimator which requires the solution of a dynamic programming problem at each time step. Estimation using the receding horizon method results in increased computational complexity but provides very little phase lag between the estimated and true excitation forces [7]. A similar disturbance model estimation method to [6] is [9], where a linear Kalman filter is used on simulation data to perform an investigation into excitation force estimation and prediction. The direct estimation method shown in [7] applies the same principle that is used for drag estimation in [10] and unknown parameter identification in [11]. The work of [12] provides an overview of control algorithms, many of which rely on the use of a linear Kalman filter. Excitation force prediction is of primary interest in the overview [12].

This work describes the implementation of an extended Kalman filter that uses a nonlinear continuous time model and discrete time measurements. The use of a nonlinear Kalman filter enables the inclusion of the nonlinear Morison drag term, that has become a commonly included in many lumped parameter models and WEC modeling software packages [13, 14]. The implementation and sensitivity analysis of wave excitation force estimation with an extended Kalman filter on experimental data from [15] is described using both the direct estimation and disturbance estimation methods. The considerations that need to be taken into account when implementing these two methods with a nonlinear model are given, and the sensitivity of the estimation methods to various conditions are investigated.

The work described in this paper is differentiated from previous literature in this field by providing the following contributions:

The excitation force estimation is performed on experimental wave tank data.

• 2.

The implementation of a nonlinear Extended Kalman Filter made use of two different methods of adapting the Kalman filter to estimate the excitation force.

• 3.

The inclusion of an estimated velocity for the water surrounding the float allows for the nonlinear Morison drag term to be used in the model.

• 4.

The sensitivity of the estimation algorithm to measurement noise, wave climate, sampling rate, radiation realization, and the number of frequencies that are used to model the excitation force are studied.

• 5.

An equal energy approach is taken to choose the wave frequencies that comprise the disturbance model of the wave excitation force.

This work provides insight into the implementation of a nonlinear Kalman filter with experimental data for wave energy conversion. Reproducible methods are described to facilitate implementation of the EKF for excitation force estimation and an investigation of the limitations of the EKF is performed to help ensure the suitability of the application. The purpose of this work is to enhance the performance of state observers when estimating the wave excitation force acting on WECs with the goal of enabling more effective feedback control.

Advanced control strategies such as latching control or Model Predictive Control (MPC) offer increased power production when compared to traditional WEC control strategies [16, 5]. Optimal power absorption is achieved when the phase of the WEC is properly aligned with the wave excitation force, to this end latching control uses an actuator to lock a WEC in a specific configuration long enough for the WEC to achieve the optimal phase with the excitation force [2]. Model predictive control and Nonlinear Model Predictive Control (NMPC) are methods that solve a series of successive, small time horizon, optimal control problems [17]. The optimal control problem is solved for a future timespan and takes into account the current wave climate when determining the optimal control input. However, a short prediction of incident wave elevation or wave excitation force is necessary in order to implement this type of controller [18]. As the wave excitation force estimation is enhanced forecasting methods can be better implemented to provide the necessary time horizons for predictive and latching control methods.

The dynamic model used in this work is described in section 2. The experimental tests used in the development and implementation of the estimation algorithms is described in section 3. A description of the identification algorithms used is given in section 4. Excitation force estimation results and sensitivity analyses are given in section 5. Finally, concluding remarks are presented in section 6.

2 Wave Energy Converter Model

The potential based time domain integro-differential equation analyzed in [19] has become the most prolific lumped parameter model for a WEC Equation of Motion (EOM). However, semi-emperical models or simplified lumped parameter models are also possible for WECs as shown in the detailed modeling description given in [10]. The equation of motion of a semi-submerged body in heave is given by

[TABLE]

where the inertial force of the object is the product of the “virtual mass”, $M_{v}$, and the heave acceleration, $\ddot{y}$. This “virtual mass” is the sum of the mass of the body and the infinite frequency added mass of the entrained fluid, as is traditionally used in the Cummins formulation [20]. The terms used in (1) are defined as: $F_{B}$ is the hyrdostatic force on the object, $F_{D}$ is the viscous drag force, $F_{R}$ is the radiation damping, $F_{M}$ is the restoring force from mooring lines, $F_{ex}$ is the wave excitation force, and the displacement of the body is $y$. The dot notation represents derivatives with respect to time, for example, $\dot{y}=\frac{dy}{dt}$ and $\ddot{y}=\frac{d^{2}y}{dt^{2}}$. Discussion of the simplification of a general equation of motion to the form shown in (1) can be found in [10].

The potential based force terms of (1) are developed in agreement with the Cummins equation [20] which was originally developed to describe the motion of ships and is now applied to semi-submerged bodies. The viscous drag term is the lone nonlinear term in (1) as described in [19, 13] and the inertial, excitation, restoring, and radiation forces are assumed to be linear for the following development.

2.1 Mooring and Hydrostatic Forces $F_{M}$ and $F_{B}$

The forces resulting from mooring and umbilical cables are modeled together as the product of a constant parameter and the body displacement as given by the linear relationship

[TABLE]

where $K_{mu}$ is lumped linear spring constant describing the forces acting on the body by the mooring lines and umbilical cable and $y$ is the displacement from the static equilibrium point. By modeling the response of the body from the equilibrium position, as with the hydrostatic restoring force, only a linear relationship is needed to determine the vertical force contribution of the cables as a function of heave position. The constant can be determined experimentally using load cells and displacement measurements, described in section 3

The hydrostatic restoring force is determined by

[TABLE]

where $A$ is the projected area that is being submerged, $\rho$ is the density,$g$ is the gravitational acceleration constant, and $y$ is the displacement from the static equilibrium. It is important to note that (3) assumes that a constant area is being submerged. While a linear force was used in this work, nonlinear buoyant forces have been considerd by [21].

2.2 Drag Force $F_{D}$

The viscous drag is modeled using a quadratic Morison drag form [22, 23]. The drag in the heave direction is given by

[TABLE]

where $V$ is the relative velocity between the object and the water particles surrounding it. The density of the fluid is given by $\rho$, and $A_{p}$ is the projected area of the object in the heave direction. A signed velocity square term is computed by taking the product of velocity, $V$, and the absolute value of the velocity, $|V|$. Finally, the drag coefficient is given by $C_{D}$. A discussion on the treatment of the viscous drag force as well as an identification method can be found in [10] A small linear damping force $bV$ is included to compensate for the fluid shear forces acting on the body and linear damping in the mooring lines. The coefficient $b$ is based on the experimental response of body at low velocities.

2.3 Radiation Force $F_{R}$

The force acting on the WEC that results from radiated waves is calculated in the time domain through the use of the convolution integral as illustrated in [20]. A state space realization of the convolution integral is created as an efficient method of computing the radiation force. The state space representation of the radiation force can be formed as

[TABLE]

where $X_{r}$ is the state variable of the realization. The state matrices $A_{r},B_{r},C_{r},$ and $D_{r}$ are a non-unique linear realization of the convolution integral and are generated from the impulse response function [24] of the system. A more detailed description of the Fourier transform inversion and computational considerations of describing the radiatoin force as a state space system is given in [10]. The work of [25] states that for the radiation realization to have a causal input-output relationship, $D_{r}=0$.

2.4 Excitation Force $F_{ex}$

The wave excitation force can be calculated using the excitation force Impulse Response Function (IRF), $F_{IRF}(t)$, and the wave surface elevation, $y_{w}$, using the convolution integral,

[TABLE]

This convolution is used to compute the true value of the excitation force for this work. The time range of the numerical integration is truncated to be longer than the active response of the excitation IRF which is determined by applying the inverse Fourier transform on the frequency dependant excitaiton force calculated from a BEM [26].

3 Available Data

Experimental data was made available by the Oregon State University Northwest National Marine Renewable Energy Center’s Marine and Hydrokinetics Open Data Project [27]. The data set that is used in this work was a physical model test of a 1:10 scale floating power system at the COAST lab at Plymouth University [15]. The deep water tests described in [15] is used to determine the parameters in the equations of motion.

The mooring study described in [15] used a wave tank to compare several mooring configurations of two bodies in water. Of these two bodies, the scaled Floating Power System (FPS) is analyzed in this current work as it provides a very similar geometry to the float of a point absorbing WEC. For this work a subset of experimental runs were chosen. Seventeen different panchromatic runs using an elastic “o-ring” as a taught mooring are used to test the Kalman filter estimation algorithm in section 5. Table 1 shows the significant wave height, $H_{s}$ and peak period, $T_{p}$ of the incident wave spectra. Each data set had roughly six-hundred seconds of data collected.

Incident wave elevation and mooring line forces were measured at 128Hz and motion data was recorded at a sampling frequency of 200Hz. Irregular waves were generated using the Bretshneider spectrum.

3.1 Reproducing Numerical Work

Detailed descriptions of the experiments along with wave tank schematics and dimensioned drawings can be found in the Marine and Hydrokinentics Open Data Project [27]. Key paramters used in the WEC models not found in [15] or [27] can be found in [10]. Frequency dependent terms, such as added mass, radiation damping, and wave excitation are calculated using the BEM code NEMOH [28]. The scripts used for this paper to execute NEMOH and all other figures will be made available on GitHub [29].

4 Continuous-Discrete Extended Kalman Filter

The Kalman filter incorporates measured outputs, a model, and Gaussian error description to determine a statistically significant estimate of the state dynamics.

A nonlinear continuous time model with discrete time experimental data points can be observed using a Continuous-Discrete Extended Kalman Filter. The EKF is a nonlinear adaptation of the Kalman filter which observes a dynamic system in which zero-mean Gaussian white noise is assumed in both the dynamics and measurement. The model that is observed by the extended Kalman filter is described by

[TABLE]

where f is the nonlinear state dynamics, h is the nonlinear output function, $\textbf{G}(t)$ is a vector that describes the noise of each state, $\textbf{w}(t)$ is zero-mean Gaussian white noise with a variance $Q(t)$, and $\textbf{v}(t)$ is zero-mean Gaussian white noise with a variance of $R(t)$. The states, $\textbf{x}(t)$, and output, $\tilde{\textbf{y}}(t)$, are written in bold font to symbolize vector quantities. The tilde symbol, $\tilde{\bullet}$, modifying the output, y, represents that this is an experimentally measured quantity. The accuracy of the Kalman filter relies on knowledge of the process and measurement noise terms. Correct quantification of the process noise is essential for the EKF to function effectively. For discussions concerning the determination of the process noise covariance see [11, 30]. Since the EKF deals with nonlinear systems that must be integrated to determine state dynamics, reasonable initial conditions must be used to ensure filter stability [11].

The Kalman filter will iterate through the gain, update, and propagation steps for each discrete time in order to estimate the states of the system. The gain step calculates a matrix of gains that is used to correct a previous step’s estimates based on the covariance of each state. The update step uses the calculated gain to adjust the covariance of each state and estimated value based on the difference between estimated and measured states. Finally the propagation step calculates the covariance and state estimates of the next time step.

Details on the theory of the continuous-discrete EKF can be found in [11] and the implementation of the algorithm on a WEC application can be found in [10]. For this implementation of the EKF the Joseph stabilized form of the covariance update equation is used, and to ensure the symmetry of the covariance matrix is enforced by taking the element-wise average of the covariance and its transpose at each step. In order to use the Kalman filter for estimation the model must be modified to provide the estimated parameter as a function of the observed states as described below.

4.1 Direct Identification Estimation

To directly identify a force or parameter using the EKF the variable in question is appended as an extra state with no dynamics. Each state in the dynamics of a Kalman undergoes a “random walk” to determine the best estimate of the state variable at each time step [11]. The nonlinear dynamics in (8), described in section 2, when reduced to a first order system is,

[TABLE]

where $x_{r}$ is a vector of radiation states, $x_{1}$ and $x_{2}$ are the first order variables for position and velocity respectively, and both $F_{ex}$ and $\dot{y}$ are inputs to the system representing wave excitation force and water velocity, respectively. The estimated wave excitation force can be used to provide a functional estimate of the water velocity surrounding the WEC as described in section 4.3.2.

The inclusion of $F_{ex}$ as an estimated force necessitates the inclusion of another state in the dynamics as,

[TABLE]

where $F_{ex}$ is now the third state of the system estimated by the EKF. The inclusion of a state with zero dynamics means that the “random walk” that takes place during the gain and propegation steps of the Kalman filter is responsible for the changes in the estimated value. Equation (11) assumes that the time derivative of the excitation force has the form $\dot{F}_{ex}=gw(t)$, where $g$ is a single term of $G(t)$ in (8) and $w(t)$ is a Gaussian zero-mean white noise signal. The EKF is then modeling the excitation force simply as noise where the current value of $F_{ex}$ does not depend on previous estimates. The choice of the process noise $w(t)$ for the estimated force is crucial to success of the estimation algorithm, the selection of noise is discussed in [11] and the choice used in this work is described in section 4.3.1.

The direct identification method only adds a single state to the dynamics of the system and as such is the most computationally simple identification method used in this work. An additional initial condition is necessary with the addition of any state is required, however precise determination of this initial condition is not necessary as it is expected that there be an initial error in the dynamics. This method of identifying the excitation force does not require information on the wave climate, however a sensitivity analysis of the process noise corresponding to the identified state should be performed. An example where a large process noise for an identified state is able to identify a time-varying viscous drag coefficient is shown in [10].

4.2 Disturbance Model Estimation

An alternative method of estimating the wave excitation force is to model the incident wave force as a disturbance to the WEC dynamics [31]. Ocean waves and excitation forces are typically modeled as harmonic oscillator, or the sum of a large number of sinusoids, with parameters determined by a wave energy spectrum and randomized phase [4]. The disturbance, as described by [6], will be modeled using a similar approach. The regular wave representation of the excitation force is given as

[TABLE]

where $\mathbb{R}$ represents the function to extract the real part of the argument, $H$ is the peak to peak wave height, and $F_{X}(\omega_{r})$ is the magnitude of the frequency domain excitation force for a given frequency, $\omega_{r}$. The imaginary unit is donated by $i$, and $t$ is the time in seconds. This excitation force can be generalized to irregular waves by

[TABLE]

where $S(\omega)$ is the power of the wave spectrum as a function of frequency. The variable $\Delta\omega$ is the discrete distance separating each frequency and $\phi$ is vector of randomly generated phase angles in radians [32].

Equation (13) shows an infinite integral over a range of frequencies which is typically modeled a discrete sum [4]. By looking at a single frequency component of the wave excitation force a description of the excitation force dynamics can be shown as,

[TABLE]

The differential equation describing the excitation force can be implemented in the first order system of equations by

[TABLE]

where $f_{ex}$ and $\dot{f}_{ex}$ are now states of the dynamic system. The above expression describes a system where a single excitation force signal is estimated, similar to the direct identification method. However, (15) enforces that the estimated excitation force be a solution to the differential equation given in (14).

This excitation force estimation method can be generalized to use multiple wave frequencies by using

[TABLE]

for the dynamics. The bar over the excitation force denotes a vector of disturbance forces. The excitation force becomes $F_{ex}=\Sigma\bar{f}_{ex}$ to denote the sum of $N$ excitation force terms that each oscillate at a different frequency. The term $\Omega$ denotes the diagonal matrix containing each of the $N$ frequencies corresponding to the entries of $\bar{f}_{ex}$. The sum of each excitation force frequency allows for the excitation force to be modeled in a way that realistically describes the time derivative behavior of the wave excitation force. This method however requires knowledge of the wave climate for an appropriate choice of the disturbance wave frequencies. The choice of these frequencies is discussed in section 4.2.1.

Just as with the direct identification method, each of the new states that are included in the dynamics are subject to a “random walk” that is influenced by the process noise covariance in the gain and propagation steps of the Kalman filter. The disturbance estimation method becomes nonlinear if any of the frequencies in $\Omega$ are varied, where any frequencies would then be appended to the dynamics as an extra state [6]. For this work the frequencies remain constant.

The initial conditions of the excitation terms should satisfy the differential equations given by (14). These initial conditions act as the terms that dictate the phase of the excitation force terms, however after a small amount of simulation time the Kalman filter is able to adjust to reasonable errors given by the initial conditions. As a result, once the Kalman filter is properly set up then nearly any physical initial condition can be used to start the estimation.

The disturbance model introduces two states for every disturbance frequency that is used in the dynamics which is a significant increase in computational complexity over the direct identification method. However, the use of the equal energy approach when selecting wave frequencies provides an effective way to model the incident wave with only few frequencies. This enables the use of the disturbance model to more realistically describes the wave excitation force.

4.2.1 Equal Energy Approach

The implementation of the disturbance method requires the choice of $N$ frequencies to comprise the wave excitation force. The typical approach to representing irregular waves uses many frequencies, on the order of 10,000 sinusoids [7], however this number would be computationally impractical for implementation in a Kalman filter. The wave frequencies used to construct an irregular wave can be found in several ways as described in [33]. Some examples include an arithmetic progression, which selects wave components that are equally spaced in frequency, and a geometric progression, which selects wave components that have a constant ratio of wave frequency. Both of these methods typically require large numbers of waves to realistically represent an irregular wave. In contrast, the equal energy method can be used to realistically represent an irregular wave with relatively few frequency components, and is mentioned in [4]. The basic premise of this method is to select wave components such that each frequency represents an equal amount of spectral energy. Since this energy is calculated by considering the area under the energy spectrum, this method is sometimes referred to as the equal area approach.

Figure 1 shows the wave energy spectra of the 17th run of the available experimental data described in section 3. The significant wave height of the generated wave was 375 mm and the peak period was 2.37 seconds. The vertical blue lines show the divisions between each equal energy area of the spectra and the red circles indicate the power and frequency of the components that were selected to form the disturbance equations.

The divisions are created by determining the area under the energy spectra and dividing this by the number of desired frequency components, which was in this case 5 components, to determine the area of each division. A threshold for minimum power is required otherwise the very high and very low frequency components would artificially bias the selection towards frequencies with very little energy. Once the equal energy areas are determined a single frequency with the average power is selected from each equal energy area to provide the average power for that division. This method of determining the disturbance frequencies is effective for selecting wave frequencies and initial estimations of each components magnitude while only requiring the peak period and significant wave height of the incident wave spectra.

4.3 Implementing the EKF

In order to implement the extended Kalman filter a few key considerations were required prior to successful estimation of the excitation force. First, and often the most difficult part of any Kalman filter, is determining the process noise covariance given in section 4.3.1. The other barrier to the implementation of the nonlinear equation of motion is that in typical WEC models both the excitation force and incident water velocity are unknown, therefore an approximation of the wave velocity is required. The method used in this work for estimating the water velocity is described in section 4.3.2.

4.3.1 Process Noise

The process noise covariance is the product of $G(t)\textbf{w}(t)$ that describes the error in the dynamics of the system. The diagonal entries of the symmetrical matrix $G(t)\textbf{w}(t)$ are variances of the noise in each state and non diagonal entries quantify the correlation between the states. A typical assumption for the development of a Kalman filter is to assume that no cross-correlation exists between states [11]. For the equation of motion the radiation states are excepted from having a process noise variance as the error in the radiation force is included in the error for the position and velocity states [10].

The process noise for each state was found by simulating the response of the model to the inputs determined from experimental measurements. The convolution formulation of the excitation force in (7), is used to determine the wave excitation force that the float experienced during each experimental run, and the water particle velocity is known from measurements of the incident wave 3. The dynamics of the float are simulated using the experimental measurements and the states are compared to the measured states from [15]. The process noise describes the error between the simulated states and the experimental states, so the variance of this simulation error is calculated for each run and used as the process noise covariance in the EKF.

Recall that we are considering various experiemental runs from [15], and attempting to determine the wave excitation force for each run. Since the process noise is a variable, with dimensions, that will vary based on the incident wave it becomes necessary for each run to have its own process noise. A typical Kalman filter implementation requires a significant amount of tuning to provide efficient estimation is achieved [11] and in some cases optimization algorithms are used to identify a “best” choice of the process noise terms [31]. For this work, a reasonable value of process noise for the excitation force was chosen for each separate sensitivity analysis as changes in sampling rate, model complexity, and measurement noise all change the robustness of the EKF to changes in the estimated state process noise.

4.3.2 Estimating Water Velocity

The estimation of the water velocity of the incident wave enables the calculation of the Morison drag force as the relative velocity is needed between the body and the surrounding water particles. With this estimation of the velocity it becomes possible to observe the model described in section 2 using an EKF. It is important to note that this nonlinear time domain model is well represented as a generic model for point absorbing WECs [13].

The estimation of the water particle velocity relies on a regular wave representation of the wave excitation force using (12). When simulating the excitation forces in irregular waves using (7) it is necessary to incorporate the time delay and the frequency dependent magnitude of the excitation force, however in order to approximate the incident wave position the magnitude is assumed to have a small phase lag and a magnitude equal to the magnitude of $F_{X}(\omega_{p})$, where $\omega_{p}$ is the peak frequency of the exciting wave.

By using the linear relationship of (12) between the wave elevation and excitation force, the wave elevation can be approximated by

[TABLE]

where $F_{ex}(t)$ is the estimated excitation force and $||F_{X}(\omega_{p})||_{2}$ is the L2-norm, Euclidean norm, of the excitation force at the peak frequency of the incident wave. This norm is required as each entry of $F_{X}(\omega)$ is an imaginary number that is used to represent the magnitude and phase of the wave excitation force. The wave velocity can then be computed using a finite difference method based on the previous estimations of the water elevation. This simple method allows for the full implementation of the Morison drag equation with the time domain model discussed in section 2. The estimation results of the water particle velocity are given in section 5.

5 Estimation Results

The two estimation methods presented in this work were tested using the experimental data described in section 3. The discrete measurements used in the EKF are the heave position and velocity, and the estimated excitation force is compared to the excitation force calculated using the non-causal convolution integral formulation in (7). This convolution integral relies on knowledge of the time history of the incident wave as well as the impulse response of the buoy in the heave direction, both of which are known before estimation as a result of the work in [15]. The experimental values of water velocity is determined by taking the time derivative of the measured water surface elevation.

Each estimation run is evaluated using the Normalized Mean Square Error (NMSE), given by the expression

[TABLE]

where $s(t)$ is the time domain reference signal, the bar notation denotes the mean value, $\tilde{s}(t)$ is the estimated signal, and the notation $||\bullet||_{2}$ represents the L2-norm of a vector. The NMSE goodness of fit parameter is a scalar with a value less than one. Alternatively the goodness of fit can be represented as a percentage as shown in [7]. In this work the value will be represented as a scalar where $NMSE=1$ would denote perfect estimation of the reference signal [34] . In the case of this work both the incident wave velocity and excitation force are evaluated using (18).

Figure 2 shows the time domain estimation of the wave excitation force using both the direct identification method and the third order disturbance estimation method, where third order signifies that the disturbance models is a sum of three excitation force frequencies, for the incident wave. The reference signal that was calculated from experimental measurements is shown in black, and the results from the direct identification and the disturbance estimation methods are shown in blue and red, respectively. The third order disturbance estimation is shown as it provides a representative example of the performance of the disturbance estimation method while also being directly comparable to the work of [6]. This figure shows that some of the very sharp changes in the excitation force are not captured by the estimation algorithm, however it is apparent that the estimation algorithm results in a excitation force signal that represents well the true signal. For this estimation run the significant wave height of the incident wave was $H_{s}=300$ mm and the peak period was $T_{p}=2.37$ s. The NMSE of direct and disturbance estimations are $0.914$ and $0.927$, respectively. In order to test the estimation algorithms in a variety of wave climates the excitation force estimation was performed for the seventeen different wave climates described in section 3, the results of which are shown in Fig. 4 which is described below.

The estimation of the incident wave water velocity is shown in Fig 3 for the run 8 in table 1. The water velocity results from the direct estimation implementation and the disturbance estimation implementation are once again shown in blue and red, respectively. The method by which the incident water velocity is estimated is described in section 4.3.2. It is not difficult to tell by visual inspection that estimation of the water velocity is not as accurate as the estimation of the excitation force. This is reflected in the goodness of fit values, which for the direct and disturbance estimation methods results, shown in Fig. 3, are $0.825$ and $0.831$, respectively. The approximation used to describe the velocity as a function of excitation force is very simplistic and a more detailed formulation of water velocity as a function of excitation force using the relationship given by (13) could be implemented. Additionally a higher order finite difference scheme could be implemented. However, as the influence of the water velocity on dissipation forces does not drive the dynamics of the WEC, this simple approximation method is sufficient for this work.

The NMSE of several estimation algorithms applied to each run in table 1 are shown in Fig. 4. Recall that both of the estimation methods being presented are run seventeen times, once for each of the different wave climates being considered. The results are shown in Fig. 4 as box and whisker plots to give an understanding of the overall performance of each method. For this comparison the direct estimation method is shown compared to the disturbance estimation method of varying orders (one excitation term to five excitation terms). It can be seen from the figure that the performance of the disturbance estimation algorithm provides the best estimation result using three excitation force components. As a result the third order disturbance estimation method will be used when analyzing the sensitivity of the estimation algorithms. The estimation of the incident wave velocity is shown in Fig. 5. As shown in this figure the mean of the seventeen NMSEs tend to be roughly 10% lower than the excitation force estimation shown in Fig. 4.

Both figures 4 and 5 show a single outlier marked with a red $+$. This outlier corresponds to the first run in table 1. This high frequency wave is shown to be consistently difficult for the EKF to estimate. While it is possible to tune the parameters of the Kalman filter to better estimate the high frequency wave that is currently being shown as an outlier, the purpose of this work is to show general trends in the estimation capabilities of the EKF methods. Therefore, rather than producing the highest possible goodness of fit for a single wave climate we have focused on tuning the parameters of the Kalman filter to produce acceptable estimation for a wide range of wave climates. It is of interest to note that this first run consisted of a wave that generated relatively small excitation forces when compared to the waves from the other runs being considered. For example, the excitation force of the seventeenth run has a magnitude about four times larger than the first run. This means that the wave from the first run will have a much smaller effect on the dynamics of the buoy than any of the other runs considered in this work. Therefore, while this first run consistently has the lowest NMSE of any of the runs being considered, the waves that are consistently estimated well by the EKF have a much larger effect on the buoy dynamics.

As each run represents a different incident wave the NMSE can be plotted against several parameters as shown in Figures 7, 6, and 8 which show the NMSE as a function of peak period, significant wave height, and size ratio respectively.

In Fig. 6 the NMSE is shown to not be directly influenced by the significant wave height, unlike the wave frequency and size ratio. The performance of the estimation algorithms tend to be worse in the cases where the wave is oscillating at a high frequency. Figure 7 shows that the lowest performing estimation runs occur at low frequencies. The size ratio is commonly used to determine the dominant forces for modeling purposes [35]. As the relative size of the float when compared to the wavelength becomes large the assumptions used in the development of the EOM for this work may be called into question as the size ratio is shown to negatively impact the NMSE of each run in Fig. 8.

5.1 Sensitivity to Sampling Rate

The motion data is originally recorded at a sampling frequency of 200Hz. A subset of the motion data is used to test the sensitivity of the estimation algorithms to sampling rate. Figures 9 and 10 show the response of the direct estimation and third order disturbance estimation when the sampling rate is lowered. To select a subset of the measurements an evenly spaced number of the measurements were selected and treated as the measured data points. The resulting sampling frequencies are used as the control variable to determine the influence on the NMSE of the estimations.

Figures 9 and 10 show that even at much lower sampling rates the effectiveness of the EKF to identify the excitation force is still resulting in a NMSE greater than 0.9. The outlier corresponding to the first, high frequency, run shows a more significant performance degradation than the mean performance of the rest of the runs. It intuitively follows that a low sampling rate would adversely affect the estimation of a high frequency wave. It is important to note that while the average for the runs appears to increase at a sample rate of 2Hz the fact that an extra outlier is shown means that the overall performance for all of the runs is lower. This extra outlier in both figures 9 and 10 corresponds to run 4 which has a low period for its amplitude. The robustness of the estimation algorithms to a change in sampling rate shows promise for the potential of real time estimation for WEC deployments.

5.2 Sensitivity to Radiation Force Order

The radiation force is modeled using a state space realization given by (6) that can vary in size. Reference [31] shows a study in which the order of the radiation realization is varied and the effect on the estimator performance is evaluated using simulated wave measurements. When generating the state space realization of the radiation force it is typical to use a coefficient of determination as shown in [10]. The fourth order state space realization, used in this work, was the smallest realization that satisfied the minimum coefficient of determination threshold established in [10]. However, each state used in the radiation realization increases the computational burden of the model. Therefore, a simple investigation into the effects that the order of the radiation force has on the excitation force estimation is conducted.

It is shown in both Figs. 11 and 12 that using a first order realization for the radiation results in a significant decrease in performance. Note that for the second through fifth order realizations the performance of the estimator remains fairly constant. However, as mentioned above, it was shown in [10] that the fourth order radiation force is the smallest realization that is viable for identifying a missing drag force for this heaving float. As a result, the fourth order radiation force is used throughout the rest of the simulations in this work.

5.3 Sensitivity to Measurement Noise

The experiments described in [15] obtained high quality position measurements by using a Qualysis data acquisition system. The measurement noise used in the EKFs is very small, calculated from the residuals of the position measurements. The standard deviation of the measurement noise for position and velocity was $v=6.075\times 10^{-4}$. To investigate the sensitivity of the estimation algorithms to measurement noise an artificial, zero mean, Gaussian noise term was added to the measurements. Figures 13 and 14 show the NMSE for the direct estimation and disturbance estimation methods, respectively, as artificial measurement noise is added to the experimental measurements. The actual measurement noise is very small, the artificial measurement noise, $v_{a}$ is described by a noise factor that is an integer multiple of the experimental measurement noise. Both Figs. 13 and 14 show the resilience of the EKF, as no appreciable performance degradation is shown until the artificial noise is several orders of magnitude larger than the measurement noise.

6 Conclusion

This paper presents the implementation of an Extended Kalman Filter (EKF) to estimate the wave excitation force of a Wave Energy Converter (WEC). Many advanced control strategies require the wave excitation force that is acting on the WEC to determine the control input. However, real time measurement of the wave excitation force is typically not viable for ocean deployments. To provide information about the wave excitation force, Kalman filters are implemented as optimal state observers which are configured to estimate the wave excitation force along with the dynamics of the WEC. Prior to this work, typical implementations of the Kalman filter have been restricted to linear Kalman filters and have been tested with simulations rather than experimental data.

Two different estimation methods, making use of the nonlinear Kalman filter, are described in this work, both of which are tested with experimental data from a 1:10 scale wave tank test of a floating power system which closely resembles a point absorbing float. The direct estimation method relies on modeling the excitation force as a single state with no dynamics, allowing for the “random walk” of the Kalman filter to determine a statistically significant estimate of the excitation force. The disturbance estimation method models the wave excitation force as a sum of sine waves with the appropriate derivatives modeled as states. An equal energy method is used to determine the amplitudes and frequencies of the harmonic components. The quantity of the harmonic excitation terms is varied and discussed in this work.

The potential based time domain integro-differential equation of motion, used in this work, is a widely applied lumped parameter model of WEC bodies. This model includes a Morison drag term which is nonlinear and requires knowledge of the incident wave velocity to describe the drag acting on the object. The nonlinearity of this model requires the use of the EKF rather than the simpler linear Kalman filter. Much like the wave excitation force, the incident wave velocity is typically unknown in WEC deployments. A simple approximation of the incident wave velocity is developed and tested in this work, which enables the full implementation of the nonlinear equation of motion. The normalized mean square error between the estimated excitation force and true excitation force is evaluated using seventeen different wave climates. The sensitivity of the excitation force estimation error is evaluated for variation in the model of radiation force, measurement noise, sampling rate as well as wave climate.

Both of the estimation methods discussed in this work were shown to provide good estimates of the excitation force using experimental data that would be available during a WEC deployment. The harmonic oscillator model of the wave excitation force utilized in the disturbance estimation method is a more realistic model of the force, which results in some improvements in estimation accuracy. However, the disturbance model adds greatly to the computational complexity of the extended Kalman filter while providing modest improvements. The excitation force estimation of both methods is shown to be quite resilient to a wide range of situations. This work shows that a common, nonlinear, lumped parameter model, previously unused for wave excitation force estimation, can be used to accurately estimate experimental excitation forces through the use of an extended Kalman filter.