1

Comparison of frequency estimation methods for
wave energy control

arXiv:1804.02247v1 [cs.SY] 6 Apr 2018

Paula B. Garcia-Rosa, Member, IEEE, John V. Ringwood, Senior Member, IEEE, Olav B. Fosso, Senior
Member, IEEE, and Marta Molinas, Member, IEEE

Abstract—A number of wave energy controllers tune the power
take-off (PTO) system to the frequency of incident waves. Since
real ocean waves are non-stationary by nature and not defined
by a single frequency component, the PTO can be either tuned
at a constant frequency characterized by the local spectrum,
or continuously tuned to a representative wave frequency. In
either case, a time-frequency representation of the waves is
expected since the wave profile changes over time. This paper
discusses about the PTO tuning problem for passive and reactive
controllers, in real waves, by comparing different methods
for time-varying frequency estimation: the extended Kalman
filter (EKF), frequency-locked loop (FLL), and Hilbert-Huang
transform (HHT). The aim is to verify the impact of such methods
on the absorbed and reactive powers, and the PTO rating. It is
shown that the mean estimated frequency of the EKF, and FLL,
converges respectively to the mean centroid frequency, and energy
frequency, of the excitation force spectrum. Moreover, the HHT
mean frequency has no correlation with the spectral statistical
properties. A comparison of the energy absorbed shows that up
to 37% more energy is obtained with the HHT over the other
estimation methods. Numerical simulations are performed with
sea elevation data from the Irish coast.
Index Terms—Energy harvesting, wave power, frequency estimation, control systems.

I. I NTRODUCTION

I

N order to optimize the energy extracted from waves, a
number of control algorithms for wave energy converters
(WECs) tune the power take-off (PTO) system to the frequency
of incident waves, see, e.g., [1]–[6]. For regular waves, the
PTO tuning is well defined since the waves are characterized
by a constant frequency sinusoidal motion. Nevertheless, real
ocean waves are non-stationary by nature and not defined by
a single frequency. Thus, the PTO can be either tuned to a
constant representative frequency of the local wave spectrum,
or continuously tuned to a time-varying frequency.
The PTO tuning consists of adjusting the damping of the
system, if a passive control method is adopted, or it may require adjustment of both the system damping and the stiffness
for a reactive control method. The tuning strategies usually
require knowledge of the incident wave frequency [1]–[5] or
the excitation force frequency [7], [8]. Since these frequencies
change with time for non-stationary signals, a time-frequency
Paula B. Garcia-Rosa and Marta Molinas are with the Dept. of Eng.
Cybernetics, Norwegian Univer. of Science and Technology, Trondheim,
Norway (emails: p.b.garcia-rosa@ieee.org; marta.molinas@ntnu.no).
John V. Ringwood is with the Centre for Ocean Energy Research, Maynooth
University, Maynooth, Ireland (e-mail:john.ringwood@nuim.ie).
Olav B. Fosso is with the Dept.of Electric Power Eng., Norwegian Univ. of
Science and Technology, Trondheim, Norway (email: olav.fosso@ntnu.no).

representation for the waves is expected – i.e. the frequency
is a function of time.
By adopting different methods to estimate the time-varying
frequency for irregular waves and real ocean waves, some
studies show that continuously tuning the PTO result in greater
energy absorption than tuning the PTO to a constant frequency
of the wave spectrum [2], [8]. In [2], two methods have been
adopted to estimate the on-line dominant wave frequency. One
is based on the sliding discrete Fourier transform, while the
other uses an analysis of the low-pass filtered incident wave,
where wave characteristics, such as the zero-up crossing and
crest-crest periods are calculated. In [8], an estimation of the
instantaneous frequency of the excitation force is obtained
by means of the Hilbert-Huang transform (HHT), and such
information is used for tuning the PTO damping of a WEC
on a wave-by-wave basis. Other methods have adopted the
extremum-seeking approach, where knowledge of the wave
frequency is not needed for tuning purposes [9], [10]. In such
cases, the PTO parameters are adapted on an hourly basis
(according to sea states variations) rather than a wave-by-wave
basis.
Furthermore, a number of control algorithms also rely on the
estimation of a time-varying frequency or the energy period of
the waves [6], [7], [11]–[14]. An on-line estimate of the excitation force frequency is obtained with the extended Kalman
filter (EKF) in [7]. In the proposed controller, the velocity
reference is set as the ratio between the excitation force and the
radiation damping, which is tuned to the estimated frequency
[7]. An adaptive vectorial control approach is proposed in [12],
where a frequency-locked loop (FLL) is adopted to estimate
the frequency of the WEC velocity.
The EKF-based method is based on a harmonic model
with one variable frequency [7]. Then, the EKF tracks only
a single dominant frequency. In contrast, the FLL method is
based on an adaptive filter structure and an integral controller,
where the selectivity of the filter and the tracking frequency
dynamics can be adjusted [15]. Moreover, the instantaneous
frequency, defined by the derivative of a phase function,
is a time-varying parameter which identifies the location of
the spectral peak of the signal as it varies with time [16].
However, the instantaneous frequency has physical meaning
only for mono-component signals, i.e., signals with a single
frequency or a narrow range of frequencies varying as a
function of time [16]. In order to calculate the instantaneous
frequency of multi-component signals, the HHT method [17]
firstly decomposes the original signal into a number of monocomponent signals through the empirical mode decomposition

2

(EMD). Such a decomposition has an adaptive basis and relies
on the local characteristics of the signal. Thus, the EMD can
extract different oscillation modes present in a wave profile.
The aim of this paper to verify the impact on the WEC
performance of adopting different frequency estimation methods for tuning purposes in real ocean waves. Three methods
are adopted for frequency estimation: EKF, FLL, and HHT.
The methods are conceptually different and vary from estimating a single dominant frequency (EKF) to estimating
the instantaneous wave-to-wave frequency of the oscillation
modes present in a wave profile (HHT).
This study considers both passive and reactive controllers.
For reactive control, the total power consists of a combination
of active and reactive powers. Therefore, it is fundamental
to individually determine these two power components to
indicate the actual absorbed power, and the power that has
to be supplied by the PTO during the conversion process.
The performance of the WEC is then measured in terms of
absorbed power, peak-to-average power ratio (PTO rating), and
ratio of average reactive power and absorbed power. Numerical
simulations with real sea elevation data from the Irish west
coast are presented.
II. M ODELING AND CONTROL OF WEC S
A. Equations of Motion
This study considers a single oscillating-body represented
as a truncated vertical cylinder constrained to move in heave.
With the assumption of linear hydrodynamic theory, and
neglecting friction and viscous forces, the motion of the
floating cylinder is described by the superposition of the wave
excitation force (fe ), radiation force (fr ), restoring force (fs )
and the force produced by the PTO mechanism (fp ):
mẍ(t) = fe (t) + fr (t) + fs (t) + fp (t) ,

(1)

where x ∈ R is the vertical position of the body, and m ∈
R+ is the body mass. The restoring force is given by fs =
−Sx, where S ∈ R+ is the buoyancy stiffness. From [18], the
radiation force is calculated as
Zt
− fr (t) = mr (∞) ẍ + hr (t − τ ) ẋ(τ ) dτ ,

(2)

where mr (∞) ∈ R+ is the infinite-frequency added mass
coefficient, defined with the asymptotic values of the added
masses at infinite frequency. The kernel of the convolution
term hr (t − τ ) is known as the fluid memory term:
2
π

Z∞
Br (ω) cos(ωt − τ ) dω ,

(3)

where Br (ω) ∈ R+ is the radiation damping coefficient, and
ω ∈ R+ is the wave frequency. From (1) and (2),

0

1
he (t) =
2π

Z∞

He (ω)eiωt dω .

(6)

−∞

he is the inverse Fourier transform of the excitation force transfer function He (ω), which has low-pass filter characteristics
for floating WECs, and ζ is the wave elevation. Notice that
(6) is non-causal, since in fact, the pressure distribution is the
cause of the force and not the incident waves [19].
B. Passive Control
A generic PTO mechanism, with a damper varying in time
(Bp ∈ R+ ), is considered for the passive control (PC). Thus,
fp (t) = −Bp (t)ẋ(t) .

(7)

For the case when Bp (t) = Bp , for any time t, and for
monochromatic waves, the maximum absorption is obtained
when [19]:
p
(8)
Bp = (Br (ω))2 + (ω(m + mr (ω)) − S/ω)2 ,
where mr (ω) ∈ R is the added mass. Equation (8) is frequency
dependent, and indicates that there is an optimal damping for
each frequency when the WEC is submitted to real ocean
waves, or irregular waves with a mixture of frequencies.
Here, the PTO damping is continuously modified, and tuned
to the excitation force frequency. From (8),
p
Bp (t) = (Br (ω̂))2 + (ω̂(m + mr (ω̂)) − S/ω̂)2 ,
(9)
where ω̂(t) is the estimated time-domain frequency of the
wave excitation force. In order to examine the impact of the
time-varying frequency on the performance of the WEC, three
methods are adopted for the frequency estimation: EKF, FLL,
and HHT.

For reactive control (RC), the PTO force consists of a
damping term and a spring term:
fp (t) = −Bp (t)ẋ(t) − Sp (t)x(t) ,

(10)

where Sp ∈ R is the stiffness coefficient. For incident regular
waves, if Sp (t) = Sp for any time t, and
Sp = ω 2 (m + mr (ω)) − S ,

(11)

then the reactive part of the total impedance
[Rp + Rr (ω)] + jω[(m + mr (ω) − (S + Sp )/ω 2 )]

0

Zt
M ẍ(t)+ hr (t − τ ) ẋ(τ ) dτ +Sx(t) = fe (t)+fp (t) ,

−∞

C. Reactive Control

0

hr (t) =

with M = [m + mr (∞)]. The excitation force, i.e., the force
due to the incident waves is given by
Z∞
fe (t) =
he (t − τ ) ζ(τ ) dτ ,
(5)

(4)

is cancelled and the velocity of the floating body is in phase
with the excitation force [19]. In such a case, the PTO damping
(8) becomes
Bp = Br (ω) ,
(12)
and the greatest wave energy absorption is obtained.

3

Following the same procedure for PC, Sp and Bp are tuned
to the excitation force frequency. From (11) and (12),

structural model from [21], (19) can be written in a recursive
non-linear state space form as

Sp (t) = ω̂ 2 (t)(m + mr (ω̂(t))) − S ,

(13)

υ[k + 1] = f (υ[k]) + $[k]

Bp (t) = Br (ω̂(t)) .

(14)

Notice that, for practical application studies, the physical
limits of the body excursion and the PTO should be taken
into account. Here, a PTO force constraint is implemented as a
saturation. However, this is only a theoretical approach, since
the body motion is also a function of the excitation force,
which is an external force that cannot be manipulated. The
implementation of saturation on the force signals of a real
WEC is not physically possible.
D. Energy and Power
The mean power and energy absorbed by the WEC over a
time range T are, respectively,
Z T
E
Bp (t)ẋ2 (t)dt .
(15)
P̄ = , and E = −
T
0
Notice that, for the reactive control, the delivered power has
two components: the absorbed power (or active power, that
is the power delivered to the damping Bp ), and the reactive
power (the power delivered to the spring Sp ) [19]. The mean
reactive power, and reactive energy, are respectively,
Z T
Er
, and Er = −
P̄r =
Sp (t)x(t)ẋ(t)dt .
(16)
T
0
Since the spring term in the PTO force (10) and the body
velocity can have opposite signs in (16), the PTO has to return
energy for some parts of the wave cycle. Then, the PTO system
should be able to implement bidirectional power flow for RC.
In this study, the performance of the WEC in terms of mean
absorbed power (P̄ ) is measured by the capture width ratio,
CWR =

P̄
,
2rPζ

(17)

where r is the cylinder radius and Pζ is the transported wave
power per unit width of the wave front. In deep water [19],
Z
ρg ∞ Sζ (ω)
Pζ =
dω ,
(18)
2 0
ω
where Sζ is the wave spectrum, ρ is the sea water density, and
g is the gravitational acceleration.
III. E STIMATION OF THE WAVE EXCITATION FORCE
FREQUENCY

A. Extended Kalman Filter
In order to estimate the frequency by means of the EKF,
we assume a harmonic model with a single time-varying
frequency component and amplitude, as proposed in [20].
Thus, fe (t) can be expressed in discrete-time (t = kTs ) as:
fe [k] = A[k] cos (ω[k]kTs + ϕ[k]) + η[k] ,

(19)

where A is the amplitude of the wave excitation force, ϕ is
the phase, and Ts is the sampling time. Following the cyclical

fe [k] = h(υ[k]) + η[k]

(20)


T
where υ∈R3 is the state vector defined as υ= ψ ψ ∗ ω ,
$ ∈ R3 and η ∈ R are zero-mean independent random processes with covariance matrices defined, respectively, as
E[$$T ] = R and E[ηη T ] = Q. Functions f (·) and h(·) are,


cos (ω[k]Ts k) sin (ω[k]Ts k) 0
f (υ[k]) = − sin (ω[k]Ts k) cos (ω[k]Ts k) 0 υ[k],
0
0
1


and h(υ[k]) = 1 0 0 υ[k], respectively [20].
The EKF obtains an estimate of the state vector υ[k] based
on observations of fe [k], and on the first-order linearization of
model (20) around the last state estimate. The EKF algorithm
is summarized in Table I, where I is the identity matrix of
order 3, Jf and Jh are the Jacobian matrices of f (.) and h(.),
denoted respectively by Jf [k] = ∇f |υ̂[k|k] and Jh [k + 1] =
∇h|υ̂[k+1|k] .
TABLE I
EKF ALGORITHM .
Prediction step:
υ̂[k + 1|k] = f (υ̂[k|k])
P [k + 1|k] = Jf [k]P [k|k]Jf [k]T + Q[k]
Innovation step:
υ̂[k + 1|k + 1] = υ̂[k + 1|k]
+K[k + 1](fe [k + 1] − h(υ̂[k + 1|k]))
K[k + 1] = P (k + 1|k)Jh [k + 1]T
(Jh [k + 1]P [k + 1|k]Jh [k + 1]T + R[k + 1])−1
P [k + 1|k + 1] = (I − K[k + 1]Jh [k + 1])P [k + 1|k]

Once an estimate of the state vector υ̂[k|k] is available from
the EKF (Table I), the amplitude and frequency of the wave
excitation force are, respectively, obtained as:
q
ÂEKF [k|k] =
ψ̂[k|k]2 + ψ̂ ∗ [k|k]2 ,
(21)
ω̂EKF [k|k]

= ω̂[k|k] .

(22)

B. Frequency-Locked Loop
An adaptive filter is implemented by means of a secondorder generalized integrator (SOGI), where the FLL estimates
the frequency of the input signal. Such a scheme is termed
SOGI-FLL, and has been proposed for grid synchronization of
power converters [15]. Figure 1 illustrates the block diagram
of the SOGI-FLL, consisting of the SOGI-QSG (SOGI quadrature signal generator) and the FLL structure. The method
was first adopted, within the wave energy control context, to
estimate the frequency components of the WEC velocity [12].
In this study, the method is adopted for estimating the wave
excitation force frequency.
The SOGI-QSG acts as an adaptive bandpass filter with
two in-quadrature output signals (ξ, ξ ∗ ), where ξ ∗ lags ξ
by 90◦ . The bandwidth of the filter is exclusively set by
the gain κ [15]. The FLL is responsible for estimating the
frequency of the input signal. Notice that the in the nonlinear

4

TABLE II
EMD ALGORITHM .
fen(t)

εω

−
+

κ

+
−

X

R

X

R

ξ ∗ (t)
X

Step 0: Set i = 1; r(t) = fe (t);
Step 1: Identify the local maxima and minima in r(t);
Step 2: Calculate the upper envelope defined by the maxima,
and the lower envelope defined by the minima;
Step 3: Calculate the mean envelope m(t);
Step 4: Set h(t) = r(t) − m(t);
Step 5: If h(t) is an IMF, go to next step. Otherwise, set r(t) = h(t)
and go back to step 1;
Step 6: Set ci (t) = h(t); r(t) = r(t) − ci (t);
Step 7: If i = N , define the IMF components as c1 (t), . . . , cN (t),
and the residue as r(t). Otherwise, set i = i + 1 and
go back to step 1.

ω̂F LL

R

−γ

ξ(t)

ωe

+ +

Fig. 1. Frequency estimation by the SOGI-FLL method.

frequency adaptation loop, εω = εξ ∗ can be interpreted as
a frequency error variable, and the parameter γ represents
the gain of the integral controller [15]. The selectivity of the
adaptive bandpass filter and the tracking frequency dynamics
are respectively determined by the tuning parameters κ and γ.
From Fig. 1, the state-space equations of the SOGI-FLL are
 

  

2
−κω̂FLL −ω̂FLL
ξ
κω̂FLL
ξ˙
=
+
fen
0
1
ν
0
ν̇
 

 
ξ
1
0
ξ
=
(23)
ξ∗
0 ω̂FLL ν
ω̂˙
= −γ(fen − ξ)ξ ∗ ω̂ ,
(24)
FLL

FLL

∗

where (ξ, ν) and (ξ, ξ ) are, respectively, the state and output
vectors of the SOGI, and fen is the normalized excitation
force. The FLL state equation1 is represented by (24).

where N is the total number of IMFs, which is defined here as
log2 Ns−1 [22], Ns is the data length, and r(t) is the residue.
The dominant IMF is identified through the comparison of
the energy of the IMF signals (Eci ) with the energy of the
excitation force signal (Efe ),
Z T
Z T
Eci =
|ci (t)|2 dt , Efe =
|fe (t)|2 dt ,
(26)
0

0

where ci (t) is the i-th IMF component. The dominant component cd (t) is the IMF with the highest Eci /Efe ratio.
Finally, the Hilbert transform (HT) is applied to cd (t) [17]:
Z ∞
cd (τ )
1
dτ ,
(27)
υd (t) = P
π
−∞ t − τ
where P indicates the Cauchy principal value. Then, the
dominant IMF is represented as an analytic signal,
zd (t) = cd (t) + jυd (t) ,

C. Hilbert-Huang Transform
In the HHT method, the wave excitation force fe (t) is firstly
decomposed into N mono-component signals (IMFs) by the
EMD. Then, the instantaneous frequency of the dominant IMF
is adopted for tuning purposes [8]. Figure 2 illustrates the
block diagram of the frequency estimation by this method.

fe (t)
EMD

c1 (t)
c2 (t)
..
.

with amplitude ÂHHT , phase φ̂HHT , and instantaneous frequency
ω̂HHT , respectively estimated as


q
υd (t)
2
2
ÂHHT (t) = cd (t) + υd (t) , φ̂HHT (t) = arctan
, (29)
cd (t)
ω̂HHT (t) = φ̇d (t) .

(30)

IV. S IMULATION RESULTS
identify
dominant
IMF

cd (t)

Hilbert
transform

ω̂HHT (t)

cn (t)

HHT method

Fig. 2. Frequency estimation by the HHT method.

The EMD identifies local maxima and minima of fe (t), and
calculates upper and lower envelopes for such extrema using
cubic splines. The mean values of the envelopes are used to
decompose the original signal into frequency components in
a sequence from the highest frequency to the lowest one. The
EMD procedure is summarized in Table II.
Then, the wave excitation force can be expressed as
fe (t) =

(28)

N
X
1
i=1

ci (t) + r(t) ,

(25)

A. Hydrodynamic parameters
The same heaving cylinder adopted in [8] is considered here.
The cylinder has a radius of r = 5 m, draught d = 4 m, mass
m = 3.2 × 105 kg and resonance frequency 1.2 rad/s. The hydrodynamic coefficients of the cylinder were computed using
the boundary element solver WAMIT [23]. The added mass,
radiation damping coefficients, and the frequency response of
the excitation force are shown in [8].
B. Real sea elevation data
The wave data was collected in 2010 from a data buoy in the
Belmullet wave energy test site, off the west coast of Ireland.
The wave data, provided by the Irish Marine Institute, consists
of wave elevation records of 30 minutes sampled at 1.28 Hz.
Six wave elevation records (referred as sea states S1-S6),
with different spectral distribution, were selected for our study.
Figure 3 illustrates the wave spectra of the sea states, and

5

Table III shows the significant wave height (Hs ), the peak
frequency (ωp ), the energy frequency (ωe ), and the mean centroid frequency (ω1 ) of the spectra. The statistical parameters
√
Hs , ωe and ω1 are respectively calculated Ras: Hs = 4 m0 ,
∞ n
ωe = m0 /m−1 , ω1 = m0 /m1 , where mn = 0 ω S(ω)dω is
the spectral moment of order n. ωp is the frequency at which
the wave spectrum is maximum.
Figure 4 shows the spectral density of the excitation force
for the selected sea states. Some of the high frequency waves
are filtered out by the transfer function He (ω), as can be noted
from Figure 3. The filtering characteristics are defined by the
shape of the floating body, so that the excitation force spectra
are characteristic of the cylinder adopted in this study.

Fig. 4. Excitation force spectra for sea states S1-S6.

by the FLL (0.93 rad/s) is close to the energy frequency of
the spectrum. The frequency estimated by the EKF is nearly
constant, while the tracking frequency dynamics of the FLL
depends mainly on the selection of the parameters κ and γ.
To ensure high frequency selectivity, and accurate direct and
quadrature √
components within frequency range 0.6 to 1 rad/s,
κ is set to 2 and γ = 0.16, as discussed in [12].
1.2

1.6

1.1

1.4

1

1.2

0.9

Fig. 3. Wave spectra of real wave data from Belmullet.

0.8
600

TABLE III
S IGNIFICANT WAVE HEIGHT Hs ( M ), PEAK FREQUENCY ωp ( RAD / S ),
ENERGY FREQUENCY ωe ( RAD / S ), AND MEAN CENTROID FREQUENCY
ω1 ( RAD / S ) OF THE SELECTED SEA STATES .

Hs
ωp
ωe
ω1

S1
1.26
0.52
0.59
0.70

S2
1.43
1.22
0.94
1.18

S3
1.18
0.52
0.80
1.13

S4
1.39
0.57
0.80
1.18

S5
1.42
0.74
0.90
1.06

S6
1.62
0.93
0.97
1.08

C. Time-frequency estimation by EKF, FLL and HHT
1) Superposition of two regular waves: In order to illustrate
how the estimated excitation force frequency differs according
to the method adopted, firstly we consider a simple incident
wave defined as the superposition of two regular waves: ζ(t) =
2 cos (2π/6t) + cos (2π/8t). The energy frequency and the
mean centroid frequency of the excitation force spectra are,
respectively, ωe,fe = 0.94 rad/s and ω1,fe = 0.96 rad/s.
Figure 5 illustrates the excitation force frequency estimated
by the EKF, FLL, and HT. The Hilbert spectrum shows that the
instantaneous frequency varies from about 0.92 to 1.4 rad/s,
with the highest energy content (highest amplitude) in the
lowest frequency. The EKF tracks a single frequency (0.96
rad/s) which represents the mean centroid frequency of the
excitation force spectrum, and the mean frequency estimated

10 5
10
7.5
5
2.5

1
650

700

0.8
600

650

700

Fig. 5. Frequency estimated (rad/s) by the EKF, FLL (left) and HT (right) for
the wave ζ(t) = 2 cos (2π/6t)+cos (2π/8t). The plot on the right represents
the Hilbert spectrum.

2) Sea Elevation Data: Table IV shows the energy and
mean centroid frequencies of the excitation force spectra, and
the mean frequencies estimated by the studied methods for
sea states S1 to S6. It can be noted that the mean values of
the frequencies estimated by the EKF converge to the mean
centroid frequency of the spectra, while the mean frequency
values of the FLL converge to the energy frequency of the
spectra. However, the mean values estimated by the HHT have
no correlation with the statistical parameters obtained from the
spectra. As has been remarked by Huang et al. [24], frequency
in the Hilbert spectrum has a different meaning from Fourier
spectral analysis. In Fourier spectral analysis, the existence of
energy at a frequency means that a component of a sine (or
a cosine) wave persisted through the entire time range of the
data, whereas in the Hilbert spectrum the wave representation
is local and the exact time of such oscillation is given [24].
The estimated frequencies for sea states S1 and S2 are
illustrated in Figure 6. From the Hilbert spectrum of the first
IMF, it can be noted that the frequency ranges 0.5 − 0.6 rad/s,
and 0.6 − 0.8 rad/s, have the highest energy content, respectively, for S1 and S2. Such frequency values coincide with

6

0.65

0.8

0.6

0.7

0.55

0.6

0.5

0.5

0.45
0

300

600

900

1200

1500

2.5
2
1.5
1
0.5
0

10 5
5
4
3
2
1

0

300

600

900

1200

1500

0

1800

1800

300

600

900

1200

1500

2.5
2
1.5
1
0.5
0

1800
10 5
4
3
2
1

0

300

600

900

1200

1500

1800

Fig. 6. Frequency estimated (rad/s) by the EKF, FLLt and HHT for (a) S1 and (b) S2. The plots in the bottom represent the Hilbert spectrum of the first
IMF.

TABLE IV
E NERGY FREQUENCY (ωe ), MEAN CENTROID FREQUENCY (ω1 ) OF THE
EXCITATION FORCE SPECTRA , AND MEAN VALUES OF THE ESTIMATED
FREQUENCIES ( ω̄EKF , ω̄FLL , ω̄HHT ). F REQUENCIES IN RAD / S .
ωe,fe
ω1,fe
ω̄EKF
ω̄FLL
ω̄HHT

S1
0.54
0.55
0.56
0.53
0.62

S2
0.66
0.73
0.75
0.63
1.01

S3
0.59
0.63
0.65
0.57
0.80

S4
0.56
0.59
0.60
0.55
0.73

S5
0.72
0.76
0.78
0.70
0.90

S6
0.82
0.87
0.88
0.80
0.97

dominant sea state frequency. Nevertheless, the PTO rating and
the maximum PTO required for the HHT frequencies are also
higher than for the EKF and FLL, especially for S2 (Fig. 8).
0.4
0.3

EKF
FLL
HHT

0.2
0.1
0
S1

the estimates from the EKF and FLL, but the HHT method
identifies the time at which such oscillations occur. Notice that
the HHT analysis returns nine IMFs, and the first IMF is the
dominant component in all studied cases, as shown in [8].

S2

S3

S5

S6

Fig. 7. CWR under PC tuned at frequencies from the EKF, FLL and HHT.

20

D. Effect of the estimated frequency on the control strategy
15

In order to limit the body excursions to 2.5 m for the studied
cases, the PTO force of the PC (7) and each term of the RC
in (10) is limited to ±500 kN.
1) Passive Control: The performance of the WEC is illustrated in Figures 7 and 8, for the cases when the PC strategy
adopts the EKF, FLL or HHT methods to estimate the excitation force frequency for sea states S1 to S6 (section IV-C2).
For all the studied cases, tuning the damping with frequency
estimates from the HHT gives greater energy capture than
tuning with EKF and FLL. The highest improvement is of
a factor of 1.27 when HHT is compared to EKF, or 1.37
when compared to FLL, for sea state S2. Moreover, the lowest
differences in the CWR is obtained for S1. Such behaviour can
be explained by the different energy spectral distributions of
both sea states. S1 is characterized by a narrowband spectrum
with a single dominant swell (low frequency waves generated
in other locations), with the energy concentrated in a narrow
band of frequencies. However, S2 is characterized by a twopeak spectrum with mixed wind-sea (high frequency waves
generated by the local wind) and swell conditions, with the
energy spread over a wider band of frequencies than S1.
In such a case, a method that calculates the wave-to-wave
frequency is more beneficial for PC than a method that gives a

S4

500
EKF
FLL
HHT

10
5
S1 S2 S3 S4 S5 S6

400
300
200
S1 S2 S3 S4 S5 S6

Fig. 8. WEC performance under PC (a) Peak-to-average power ratio (b)
Maximum PTO force required.

2) Reactive Control: Figures 9 and 10 illustrate the performance of the WEC, for the cases when the reactive control
strategy adopts the EKF, FLL or HHT methods to estimate
the wave excitation force frequency. For the unconstrained
cases, the CWR is very large, in most situations, and the
HHT obtains an energy improvement of up to 2.64 over
the EKF, or 1.67 over the FLL. However, the body motion
ranges from −10 to 10 m, which is practically impossible for
a WEC with a draught of 4 m. Moreover, the PTO rating for
the HHT approaches a factor of 50 and the ratio of average
reactive power and absorbed power is almost 60% in some
cases (Fig. 10.a). Such values would require oversized PTO
equipment, which would not be a rational economic choice.
For the constrained cases, the improvement obtained from
HHT over EKF or FLL is lower, but represents a more

7

4
EKF
FLL
HHT

3
2
1

300

300

150

150

0

0

-150

-150

-300
1450

0
S1

S2

S3

S4

S5

S6

0.25
0
S1

S2

S3

S4

S5

1.25

1.25

0

0

-1.25

-1.25

EKF
FLL
HHT

50

1550

-2.5
1450
600

300

300

0

0

-300

-300

-600
1450

1500

1550

-600
1450

1500

1550

1500

1550

1500

1550

Fig. 11. Time-series of the excitation force, position, and PTO force for the
constrained RC tuned with the EKF (dashed blue line), the FLL (solid red
line) and the HHT (dashed dotted black line) (a) S1; (b) S2.

12

65

1500

600

S6

Fig. 9. CWR under RC tuned at frequencies from the EKF, FLL and HHT
(a) Unconstrained case (b) Constrained case.

-300
1450
2.5

-2.5
1450

0.5

1550

2.5

1
0.75

1500

10

35
8

20
5
S1

S2

S3

S4

S5

S6

0.6

6
S1

10 7

5

S2

S3

S4

S5

S6

0.15

2.5

2.5

0

0
0

0.4

0.1

0
0.1

0.05

1.1

2.1

3.1

0
0.1

600
10

5
0.2

1200

1800

0.5

0.7

realistic scenario: an average energy improvement of 1.13 is
obtained when the HHT is compared to the EKF, or 1.23 when
compared to the FLL. Although the PTO rating varies from
a factor of 7 to 10 for all frequency estimation methods, the
ratio of average reactive power and absorbed power are much
higher for the HHT, reaching almost 14% (Fig. 10.b) for the
studied cases.
In order to illustrate the effect of the frequency estimates on
the variables of the system, for the constrained RC, Figure 11
shows samples of time-series simulation, and Figure 12 shows
the absorbed and the reactive energy over a 30-min simulation
interval for sea states S1-S2. It can be noted that, for sea
state S2, the reactive energy required for the RC, tuned with
HHT frequency estimates is much higher than the FLL or EKF
cases. In such a case, the HHT reactive power represents about
4.2% of the total power, whereas the FLL reactive power is
about 0.7%. Nevertheless, the absorbed power is 18% greater
with the HHT than with the FLL.

10

5

600

1200

1800

600

1200

1800

6

0

-5

-5

-10

-10
0

Fig. 10. WEC performance under RC: Peak-to-average power ratio (top) and
ratio of average reactive power and absorbed power as a function of CWR
(bottom) (a) Unconstrained case (b) Constrained case.

0

6

0

0.3

10 7

5

600

1200

1800

0

Fig. 12. Energy absorbed (top) and reactive energy (bottom) over a 30-min
simulation for the constrained RC tuned with the EKF (dashed blue line), the
FLL (solid red line) and the HHT (dashed dotted black line) (a)S1; (b) S2.

E. Discussion
In the EKF, a sinusoidal extrapolation method is used to
model the excitation force as a monochromatic harmonic
process, with varying amplitude and frequency. Thus, the
EKF follows a single dominant frequency. Simulation results
have shown that the estimated EKF frequency tracks the
mean centroid frequency of the excitation force spectrum.
The performance of the SOGI-FLL depends mainly on the
appropriate selection of two design parameters: κ and γ. Such
parameters define the bandwidth of the adaptive filter and the
FLL tracking frequency dynamics [15]. For the FLL, it has
been shown that the mean estimated frequency converges to
the energy frequency of the excitation force spectrum.
The HHT calculates the instantaneous frequency of the
excitation force by decomposing the signal into a number
of IMF components. Here, we have chosen only the IMF
component with the highest energy content. This component
has a wider bandwidth than the estimates from the EKF or

8

FLL, and the CWR of the WEC is greater when the HHT
frequency estimates are used for PC and RC. Even with a low
energy content, the resonance frequency of the WEC is within
the HHT estimates, which can also explain why tuning the
controllers with these estimates result in the greatest energy
absorption.
As expected, the PTO rating for PC is higher with the HHT
than with the EKF, or FLL, in most of the cases. However, for
the constrained RC, the greatest PTO rating does not indicate
the greatest energy absorption. The average value of the PTO
rating is about the same for all methods. Still, the greatest
energy absorption, obtained with the HHT, also requires the
greatest amount of reactive energy.
Here we have assumed the wave excitation force is known
completely over the simulation interval. Both the EKF and
FLL methods give an online estimate of the frequency, provided that an estimate of the excitation force is available. In
the HHT method, we have adopted an off-line EMD algorithm,
but a few implementation studies on the HHT have proposed
real-time EMD algorithms, see, e.g., [25].
V. C ONCLUSION
This paper has shown how different frequency estimation
methods, used for controller tuning purposes, impact the
energy absorbed by the WEC, the PTO rating and the required
reactive power during the conversion process. The effect of the
estimation methods on the WEC performance depend on the
control strategy employed, the PTO system constraints, and
the local wave spectrum.
For a control strategy that relies on the information of a
dominant frequency component, such as the schemes in [7]
and [13], the EKF or FLL should be adopted. The mean
frequency estimated by such methods converges, respectively,
to the mean centroid frequency or the energy frequency of
the spectrum. Moreover, if the sea state is characterized by a
narrowband spectrum, the benefit of adopting a method that
estimates the wave-to-wave frequency is relatively small.
By adopting a method that estimates the instantaneous
frequency of the excitation force (the HHT), an average
improvement in the energy absorbed of about 18% is obtained
over the EKF and FLL methods, for the constrained reactive
control strategy. For passive control, an average improvement
of 16% is also obtained for the HHT. The greatest improvements of the HHT over the other methods are obtained for
wideband spectra. In contrast to the EKF and FLL methods,
where the bandwidth is narrow and the frequency estimates
oscillate around a dominant frequency component, the HHT
frequency estimates cover a wider range, and the location of
the dominant frequency component is identified.
In this study, the first IMF component is adopted for
the HHT approach. The frequency bandwidth of this IMF
component could be narrowed by applying techniques that
deal with mode mixing in EMD. Such an approach will be
explored in future studies.
R EFERENCES
[1] G. Duclos, A. Babarit, and A. Clément, “Optimizing the power take off
of a wave energy converter with regard to the wave climate,” Journal of
Offshore Mechanics and Arctic Eng., vol. 128, no. 1, pp. 56–64, 2006.

[2] H. Yavuz, T. J. Stallard, A. P. McCabe, and G. A. Aggidis, “Time series
analysis-based adaptive tuning techniques for a heaving wave energy
converter in irregular seas,” Proc. of the Inst. of Mech. Engineers, Part
A: Journal of Power and Energy, vol. 221, no. 1, pp. 77–90, 2007.
[3] U. A. Korde and R. C. Ertekin, “Wave energy conversion by controlled
floating and submerged cylindrical buoys,” Journal of Ocean Engineering and Marine Energy, vol. 1, no. 3, pp. 255–272, 2015.
[4] C. J. Cargo, A. J. Hillis, and A. R. Plummer, “Strategies for active
tuning of wave energy converter hydraulic power take-off mechanisms,”
Renewable Energy, vol. 94, pp. 32–47, 2016.
[5] H. Mendonça and S. Martinez, “A resistance emulation approach to
optimize the wave energy harvesting for a direct drive point absorber,”
IEEE Trans. on Sustainable Energy, vol. 7, no. 1, pp. 3–11, 2016.
[6] E. Anderlini, D. I. M. Forehand, E. Bannon, and M. Abusara, “Reactive
control of a wave energy converter using artificial neural networks,”
Inter. Journal of Marine Energy, vol. 19, pp. 207–220, 2017.
[7] F. Fusco and J. V. Ringwood, “A simple and effective real-time controller
for wave energy converters,” IEEE Trans. on Sustainable Energy, vol. 4,
no. 1, pp. 21–30, 2013.
[8] P. B. Garcia-Rosa, G. Kulia, J. V. Ringwood, and M. Molinas, “Realtime passive control of wave energy converters using the Hilbert-Huang
transform,” in IFAC-PapersOnLine, vol. 50, no. 1 (Proc. of the 20th
IFAC World Congress), Toulouse, France, 2017, pp. 14 705–14 710.
[9] J. Hals, J. Falnes, and T. Moan, “A comparison of selected strategies for
adaptive control of wave energy converters,” J. of Offshore Mech. and
Arct. Eng., vol. 133, no. 3, pp. 031 101–031 113, 2011.
[10] P. B. Garcia-Rosa, F. Lizarralde, and S. F. Estefen, “Optimization of
the wave energy absorption in oscillating-body systems using extremum
seeking approach,” in Proc. of the IEEE American Control Conf. (ACC),
Montréal, Canada, 2012, pp. 1011–1016.
[11] F. Fusco and J. V. Ringwood, “Hierarchical robust control of oscillating
wave energy converters with uncertain dynamics,” IEEE Trans. on
Sustainable Energy, vol. 5, no. 3, pp. 958–966, 2015.
[12] A. M. Cantarellas, D. Remon, and P. Rodrı́guez, “Adaptive vector control
of wave energy converters,” IEEE Trans. on Industrial Electronics,
vol. 53, no. 3, pp. 2382–2391, 2017.
[13] E. Tedeschi and M. Molinas, “Tunable control strategy for wave energy
converters with limited power take-off rating,” IEEE Trans. on Industrial
Electronics, vol. 59, no. 10, pp. 3838–3846, 2012.
[14] A. H. Sakr, Y. H. Anis, and S. M. Metwalli, “System frequency tuning
for heaving buoy wave energy converters,” in Proc. of the IEEE Inter.
Conf. on Advanced Intelligent Mechatronics (AIM), Busan, Korea, 2015,
pp. 1367–1372.
[15] P. Rodrı́guez, A. Luna, R. S. Muñoz-Aguilar, I. Etxeberria-Otadui,
R. Teodorescu, and F. Blaabjerg, “A stationary reference frame grid
synchronization system for three-phase grid-connected power converters
under adverse grid conditions,” IEEE Trans. on Power Electronics,
vol. 27, no. 1, pp. 99–112, 2011.
[16] B. Boashash, “Estimating and interpreting the instantaneous frequency
of a signal - Part 1: Fundamentals.” Proc. of the IEEE, vol. 80, no. 4,
pp. 520–538, 1992.
[17] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition
and the Hilbert spectrum for nonlinear and non-stationary time series
analysis,” Proc. Royal Society London, vol. 454, pp. 903–995, 1998.
[18] W. E. Cummins, “The impulse response function and ship motions,”
Schiffstechnik, vol. 47, no. 9, pp. 101–109, 1962.
[19] J. Falnes, Ocean Waves and Oscillating Systems: Linear Interaction
including Wave-Energy Extraction. USA: Cambridge University Press,
2002.
[20] F. Fusco and J. V. Ringwood, “Short-term wave forecasting for real-time
control of wave energy converters,” IEEE Trans. on Sustainable Energy,
vol. 1, no. 2, pp. 99–106, 2010.
[21] A. C. Harvey, Forecasting, Structural Time Series Models and the
Kalman Filter. Cambridge Univ. Press, 1989.
[22] Z. Wu and N. E. Huang, “A study of the characteristics of white noise
using the empirical mode decomposition method,” Proc. Royal Society
London, vol. 460, pp. 1597–1611, 2004.
[23] WAMIT, Inc., WAMIT User Manual Versions 6.4, 6.4PC and 6.3S, 6.3SPC, USA, 1998-2006.
[24] N. E. Huang, Z. Shen, and S. R. Lon, “A new view of nonlinear water
waves: The Hilbert spectrum,” Annu. Rev. Fluid Mech., vol. 31, pp.
417–457, 1999.
[25] Y.-Y. Hong and Y.-Q. Bao, “FPGA implementation for real-time empirical mode decomposition,” IEEE Trans. on Instrumentation and
measurement, vol. 61, no. 12, pp. 3175–3184, 2012.

