IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011
627
Wave Prediction and Robust Control of Heaving
Wave Energy Devices for Irregular Waves
Marco P. Schoen,
Senior Member, IEEE,
Jørgen Hals, and Torgeir Moan
Abstract—This
paper presents a comparison of different existing
and proposed wave prediction models applicable to control wave
energy converters (WECs) in irregular waves. The objective of the
control is to increase the energy conversion. The power absorbed
by a WEC is depending on the implemented control strategy. Un-
certainties in the physical description of the system as well as in the
input from irregular waves provide challenges to the control algo-
rithm. We present a hybrid robust fuzzy logic (FL) controller that
addresses the uncertainty in the model and the short-term tuning
of the converter. The proposed robust controller uses the energy
from the power take off as input, while the FL controller portion
is used for short-term tuning. The FL controller action is based on
a prediction of the incoming wave. Simulation results indicate that
the proposed hybrid control yields improved energy production,
while being robust to modeling errors.
Index Terms—Control
systems, energy conversion, intelligent
systems, prediction methods.
I. I
NTRODUCTION
T
HE energy absorption of wave energy converters (WECs)
is influenced by the use of the knowledge of future wave
data in the control loop. The prediction of future wave data is not
a trivial undertaking, since real waves are irregular in nature. A
number of contributions have been made in forecasting irregular
waves. For example, Forsberg [1] introduced an autoregressive
(AR) with moving average (ARMA) model with fixed time hori-
zon prediction, utilizing the maximum likelihood estimation to
describe ocean waves off the coast of Sweden. The algorithm
included the on-line recursive form with an adaptation constant
for slowly time-varying sea condition. Results of the proposed
algorithm included prediction of waves for long time horizon,
i.e., 30 s using a 0.5 s sample time, yielding predictions with
adequate results. Noise considerations were given only to the
extent that outliers are removed through the standard outlier
detection rule. Forsberg’s model requires an extensive model
order in order to capture the nature of irregular waves. Halliday
et al.
[2] looked at using the fast fourier transform (FFT) to
generate prediction of waves over distance. The distance factor
Manuscript received May 29, 2010; revised September 1, 2010; accepted
December 1, 2010. Date of publication February 4, 2011; date of current ver-
sion May 18, 2011. Paper no. TEC-00239-2010.
M. P. Schoen is with Idaho State University, Pocatello, ID 83209 USA
(e-mail: schomarc@ isu.edu).
J. Hals and T. Moan are with the Center for Ships and Ocean Structures
(CeSOS), Norwegian University of Science and Technology (NTNU), 7491
Trondheim, Norway (e-mail: jorgen.hals@ntnu.no; Torgeir.Moan@ntnu.no).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TEC.2010.2101075
is used to avoid perturbations of the buoy to the measurements.
This approach necessitates the use of directional data. The sim-
ulation results presented indicate reasonable success for long
distances. One of the drawbacks of using the approach of Halli-
day
et al.
is that the spectral estimates of the FFT-based models
are less sensitive to sharp peaks and dips compared with AR
models, [3]. In addition, the FFT is for strictly stationary sys-
tems, and cannot handle the temporal fluctuations, which is an
essential component in the characteristics of irregular waves.
WECs operating in irregular waves can be categorized into
passive- and active-tuning-based control. Relevant study has
been reported in the literature since the mid 1970s [4]. Tuning
can be achieved by adapting the WECs power take off’s (PTO)
stiffness and damping to match the incident wave frequency and
the radiation wave damping coefficient, respectively. Extensive
review of this is given by Falnes [5], which includes the discus-
sion of optimal control and phase control for WECs. Latching
control attempts to achieve the same effect by halting the de-
vice’s motion in order to align the phase of the system with
the wave motion. Study on latching control has been reported
by Korde [6], [7]. Optimal latching control is discussed in [8],
where two analytical solutions are presented: one based on the
equation of motion and the other based on optimal command
theory. Using the future wave excitation, Babarit
et al.
in [9]
investigated three strategies for latching control. Though the fu-
ture wave excitation was not computed (assumed to be known),
the study has value in terms that it showed the improvement of
the energy absorbed using the wave excitation for the latching
time computation. Optimality, in terms of energy, production of
WECs has been addressed in many instances in aforementioned
literature and elsewhere. Intelligent systems, such as genetic
algorithms (GAs) may present an avenue to avoid the local min-
ima problem, and is to some limited extend explored in this
paper.
In this paper, we will propose a hybrid Kautz/AR predictive
model as well as a pure predictive Kautz model and compare
them to a predictive AR and ARMA algorithm. The outcome
of our investigation into predictive models will be utilized to
design predictive fuzzy logic (FL)/robust controllers. The pre-
diction models are based on local wave height measurements and
assume an unperturbed wave. We also attempt to address some
of the robustness aspects in designing a controller for WECs.
The robustness issues are primarily focused on the uncertainty
of parameters describing the dynamics of the buoy and not on
improved performance. The performance issues are addressed
by a proposed simple FL controller, which is responsible for
the short-term control of the device in order to maximize the
energy production. The proposed controllers in this paper do
0885-8969/$26.00 © 2011 IEEE
628
IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011
not employ latching control technology, rather use the resistive
forces of the PTO to improve the performance of the WEC. The
paper is organized in two parts. First, we describe the pursuit
of finding predictive models that work well in a stochastic envi-
ronment are of low model order (parametric models), and have
a sufficient accuracy that allows the algorithm to be used for
control of WECs operating in irregular waves. In the second
part, we propose a robust controller and combine this with a FL
controller.
II. W
AVE
P
REDICTION
M
ODELS
A. AR Predictive Model
Consider a simple discrete time finite impulse response (FIR)
AR model to characterize irregular waves
p
output vector. Defining
X
ss
=
Φ
T
Y
AR
, the updated parame-
AR
ˆ
AR
(k + 1|k) =
P
(k + 1)
×
X
ss
(k + 1).
Intro-
ter estimate is
θ
ducing the modal matrix
Q
composed of the eigenvectors of
P
= (Φ
T
Φ
AR
)
and diagonal eigenvalue matrix
Λ
of the cor-
AR
relation matrix
P,
the transformed filter weight vector results
into
ˆ
ˆ
θ
AR
(k + 1|k) =
Q
T
θ
AR
(k + 1|k).
ˆ
Transforming
X
ss
:
P
θ
AR
(k + 1|k) =
X
ss
(k + 1)
ˆ
ˆ
X
SS
(k + 1) =
QΛQ
T
×
Q
θ
AR
(k + 1|k) =
QΛ
θ
AR
(k + 1|k)
Defining
X
ss
=
Φ
T
X
ss
, which premultiplied with
Q
T
AR
yields
ˆ
X
ss
(k + 1) =
Λ
θ
AR
(k + 1|k).
(6)
Using
x
i
as the elements in
x
ss
and
λ
i
as the corresponding
eigenvalues found on the diagonal of
Λ,
each individual trans-
formed filter weight contribution is given by
a
i
(0)
(5)
y
(k|k
−
1) =
ˆ
j
=1
a
j
y(k
−
j)
(1)
where the output
y
∈
R
n
o
x1
is the wave height with
n
o
= 1
and
is measured up to time
t
=
k
−
1,
k
is the discrete time index,
p
is the model order,
a
j
∈
R
n
o
×n
o
are the influence coefficient
matrices, and
ˆ
denotes the predicted value. One of the main
problems in fitting such a time series to the gathered data is the
determination of the model order. Based on the Wold decompo-
sition theorem, an AR model can fit a variety of data character-
istics, if a sufficiently large model order is selected [10].
To construct a variable horizon prediction formulation on the
basis of (1) the following
i-step-ahead
predictor can be estab-
lished (see [11])
p
=
x
i
,
λ
i
for 1
≤
i
≤
p.
(7)
B. ARMA Predictive Model
Adding to the AR formulation a MA part enhances the flex-
ibility of the AR model [12]. Consider a typical ARMA model
in the output error formulation
p
p
ˆ
y(k|k
−
1) =
j
=1
a
j
y(k
−
j)
+
m
=0
b
m
ε(k
−
m)
(8)
ˆ
y(k
+
i|k
−
1) =
j
=1
(i)
a
j
y(k
−
j)
(2)
where we made use of the residual, which is defined by the
ˆ
difference of the estimate and the true output, i.e.,
ε
=
y
−
y.
The corresponding predictive coefficient matrices are:
a
j
=
a
1
(i)
(i−1)
where
a
j
and
b
j
are the ARMA filter matrices. The MA asso-
ciated with the
b
j
matrices is nonlinear due to the prediction
error formulation. The model can be constructed through back-
substitution of the output into the original filter model formula-
tion. An
i-step
ARMA predictor thus can be given as
p
p
(i)
a
j
y(k
j
=1
a
j
+
a
j
+1
.
(i−1)
(3)
ˆ
y(k
+
i|k
−
1) =
−
j)
+
m
=1
b
(i)
ε(k
−
m)
(9)
m
One can easily argue that by any standard, the predictive model
introduced above possesses a rather large model order. Hence,
the corresponding correlation matrix utilized in the estimation
algorithm induces a large computational burden during the op-
eration of the filter, in particular when repeated estimations are
needed to update to the current wave conditions. This is true to
some degree also for the recursive least squares (RLSs) method,
even though no inversion of large matrices is required; the model
order dictates the number of computations needed. This prop-
erty of large model order is a general characteristic of FIR filters
due to the missing feedback term in the filter structure. We can
investigate the contribution of each filter weight to the predicted
output by computing the transformed filter weights as follows:
starting with the formulation for the estimated filter parameter
matrices (least squares)
ˆ
θ
AR
= (Φ
T
Φ
AR
)
−1
Φ
T
Y
AR
AR
AR
(4)
ˆ
where
θ
AR
is a vector containing the estimated parameter matri-
ces,
Φ
AR
is the information matrix, and
Y
AR
is the respective
where the predictive parameters for the past residual terms are
defined as
a
j
=
a
1
(i)
(i−1) (0)
a
j
+
a
j
+1
(i−1)
for
j
= 1, 2, 3,
. . .
and
a
k
= 0,
for
k
>
p.
The predictive parameters for the MA terms are defined as
b
(i)
=
a
1
m
(i−1)
b
(0)
+
b
m
+1
m
(i−1)
for
m
= 1, 2, 3,
. . .
and
b
k
= 0,
for
k
>
p.
The estimation of these parameters can be done, for example,
by following the Hannan–Rissanen algorithm [13].
C. Orthogonal Basis Function (OBF) Prediction Model
One of the main reasons AR and ARMA models require large
model orders is that they are of the type FIR filter, and they do
not have a structure that allows for inclusion of prior knowledge
into the filter before the estimation occurs. AR models belong
SCHOEN
et al.:
WAVE PREDICTION AND ROBUST CONTROL OF HEAVING WAVE ENERGY DEVICES FOR IRREGULAR WAVES
629
to the class of orthogonal basis functions (OBFs), which are the
generalizations of FIR models. The basis function of FIR models
are delayed Dirac impulses expressed with the
q-operator
(delay
operator). In the pursuit of reducing the required model order
while still describing the future waveform adequately, we pro-
pose to choose orthonormal filters
L
i
(q) such that (1) becomes
p
y(k
+
i|k
−
1) =
i=1
˜
i
L
i
(q)y(k).
c
(10)
To keep the selected OBF model linear in the parameters,
L
i
(q)
are chosen before the estimation of the filter coefficient matrices
˜
i
occurs.
L
i
(q)’s
are chosen with the purpose of incorporating
c
prior knowledge of the system to be modeled. Since the waves
exhibit an underdamped dynamical characteristic, we choose
Kautz filters for the
L
i
(q)’s.
Kautz filters allow us to incorporate
conjugate complex pole pairs into the model. Considering a
given conjugate complex pole pair
(α
±
iβ),
the Kautz filter
coefficients are computed by
c
=
2α
1 +
α
2
+
β
2
and
d
=
−(α
2
+
β
2
).
(11)
Fig. 1.
Influence of Kautz parameters on prediction.
With this, the filter
L
i
(q)
can now be computed using the
following orthonormalization process:
L
1
(q,
c, d)
=
(1
−
c
2
)(1
−
d
2
)
q
2
+
c(d
−
1)q
−
d
(12)
(13)
(1
−
d
2
)(q
−
c)
L
2
(q,
c, d)
=
2
q
+
c(d
−
1)q
−
d
Step 2:
Perform a spectral analysis and identify the frequency
ω
1
with the largest contribution or magnitude in the spectrum.
Step 3:
Compute the Kautz filter coefficients
c
and
d
based
on the identified frequency
ω
1
; also assume a damping ratio
ξ.
The computation of the filter coefficients
c
and
d
require the
real and imaginary parts of the dominant root location, which
are computed as follows:
α
=
ξω
1
and
β
= sin(cos
−1
ξ)ω
1
.
Step 4:
Estimate the model parameter matrices
˜
i
using
c
−1
T
ˆ
(i)
Θ
Kautz
=
{Φ
T
Kautz
Φ
Kautz
}
Φ
Kautz
Y
Kautz
(M
−
i)
and the higher order filter parameters by using the recursive
formula as follows:
−dq
2
+
c(d
−
1)q + 1
L
2i−1
(q,
c, d)
=
L
2(i−1)−1
(q,
c, d)
2
q
+
c(d
−
1)q
−
d
−dq
2
+
c(d
−
1)q + 1
.
L
2i
(q,
c, d)
=
L
2(i−1)
(q,
c, d)
2
q
+
c(d
−
1)q
−
d
The limits on the coefficients
c
and
d
are given as
−1
< c <
1
and
−
1
< d <
1.
where
Φ
Kautz
is the information matrix,
Φ
Kautz
(q,
c, d)
as
shown at the bottom of the page, where
M
is the data length, and
Θ
Kautz
= [˜
1
, c
2
, c
3
, . . . , c
p
]
T
Y
Kautz
c
˜ ˜
˜
= [y(p +
i
+ 1),
. . . , y(M
+
i)]
T
.
Step 5:
Compute the future output using (10).
The aforementioned procedure for the computation of the
hybrid Kautz/AR filter involves two coefficients that need to
be selected. In the following, we clarify the reasoning for exe-
cuting Step 2. For this, we use a set of simulations where the
average correlation coefficient is computed between the pre-
dicted wave height and the actual wave height. The simulations
are carried out with a hybrid model that contains four Kautz
filter terms and one AR term, i.e., the model order is five, where
p
=
5 corresponds to
L
5
=
q.
The additive measurement noise
is characterized by a SNR of 1000 (0.1% noise variance), and
waves are computed based on the Pierson Moskowitz wave
spectrum [14], with a corresponding wind speed of 55.5 km/h
[30 knots, moderate Gale]. Fig. 1(a) gives the average and stan-
dard deviation (STD) for a varying
c
coefficient, corresponding
...
...
..
.
...
L
p
(q,
c, d)y(p)
L
p
(q,
c, d)y(p
+ 1)
.
.
.
L
p
(q,
c, d)y(M
−
1)
⎤
⎥
⎥
⎥
⎦
In addition to having a purely predictive Kautz filter, one can
incorporate as well some simple AR terms, resulting into hybrid
Kautz/AR filters. The step-by-step procedure for constructing
and estimating a proposed predictive Kautz or hybrid Kautz/AR
filter is given in the following. A recursive formulation can be
developed using the RLS approach, which is left for future work.
Step 1:
Collect a sufficient amount of wave height mea-
surements
y(k),
with a fixed sampling time. Normalize the
wave height time series using
y(k)
=
y(k)/y
m ax
, where
y
m ax
=
sup{y(k)}.
⎡
⎢
⎢
Φ
Kautz
(q,
c, d)
=
⎢
⎣
L
1
(q,
c, d)y(p)
L
1
(q,
c, d)y(p
+ 1)
.
.
.
L
1
(q,
c, d)y(M
−
1)
L
3
(q,
c, d)y(p)
L
3
(q,
c, d)y(p
+ 1)
.
.
.
L
3
(q,
c, d)y(M
−
1)
630
IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011
TABLE I
C
ORRELATION FOR
D
IFFERENT
P
REDICTIVE
F
ILTERS
U
NDER
V
ARYING
N
OISE
E
NVIRONMENT
to the frequency of the Kautz filter. The statistics is based on a
set of 50 simulation runs. First, we observe that when
c
=
0, the
only portions of the filter that are nonzero are corresponding to
the AR coefficients; hence, the first entry can be considered as a
simple AR filter. The correlation of the AR filter is respectable,
but using
c
= 0,
the simulations indicate that
c
is to be selected
around 0.9. Using the step-by-step procedure outlined earlier,
the average
c
value that would be selected by determining the
dominant frequency in the wave spectrum is 0.9979. The re-
lationship of
c
to the Kautz filter is that it sets its dominant
frequency to the dominant frequency of the waves.
Considering the second parameter in the proposed wave pre-
diction model, the damping, we are interested in observing the
influence of this parameter on the predicted time series. Uti-
lizing the same simulation parameters as listed earlier, with a
Pierson–Moskowitz wave spectra corresponding to a wind speed
of 46.3 km/h [25 knots, strong breeze], the statistics provided
in Fig. 1(b) is generated. The parameter
c
is computed based on
the dominating wave frequency, which resulted in
c
=
0.9973.
Fig. 1(b) indicates that the damping term
d
has little sensitivity
to the quality of the prediction. Intuitively, we would imagine
that smaller damping fits the nature of the waves. This percep-
tion is validated ever so slightly by Fig. 1(b). No damping results
into the more accurate wave prediction results.
D. Comparative Simulation-Based Analysis of AR, ARMA,
Kautz/AR, and Kautz Predictive Filter Performance
Simulations with four predictive filters are carried out in order
to compare their performance. Table I was compiled using a sea
state of moderate gale, i.e., 55.5 km/h wind speed, a moving
prediction horizon of 20 data points corresponding to a 2-s
ahead prediction at every time instant, and a model order of 35
for the Kautz/AR hybrid filter, the AR, and the ARMA filter.
The proposed predictive Kautz filter is based on a model order
of 25 odd terms. Table I details the performance for different
additive measurement noise conditions (N) in% noise variance
and is based on 51 simulation runs.
From the table, we recognize that for the noise-free case,
all filters perform less well than when some small amount of
noise is included. This behavior can be traced back to the used
estimation algorithm.
Comparing the performance of the hybrid Kautz/AR filter
with the pure AR filter for the cases where noise is present,
one recognizes that not much improvement has been achieved
by including the Kautz elements to the AR filter. The compar-
ison with the Kautz filter is to be made with caution, since it
has fewer filter weights (25) compared to the other filters (35
filter weights), which corresponds to the objective to reduce
the number of parameters used in a prediction model, while
maintaining accuracy. The employed ARMA filter has a similar
performance as the proposed Kautz filter for the cases where
noise is included. The standard deviation (STD) and average
(Ave) correlation coefficient reported in Table I indicate that
the AR and Kautz/AR filter seem to be less receptive to noise
than the ARMA and Kautz filters. Another characteristic noted
in the simulations is that the Kautz filter, and to some degree,
the ARMA filter offer a smoother prediction curve than the
AR-based filters. Correspondingly, from observing simulation
runs, the Kautz filter has a tendency to miss smaller peaks in
the forecast. A relatively simple analysis on the effectiveness to
include additional information into the filter in order to reduce
the number of filter weights is to look at the transformed filter
weights, see (7). Fig. 2 depicts the transformed filter weight
distribution of four filters. Comparing all discussed filters, the
proposed predictive Kautz filter requires fewer terms to describe
the prediction accurately.
III. B
UOY
M
ODEL
The chosen WEC model consists of a spherical heaving buoy
reacting against a fixed reference. The model includes only the
heave mode with vertical excursion
y,
while the water depth is
assumed to be infinite, and the usual linear approximation of
the hydrodynamic forces are used. Assuming that the gravity
force is balanced by the average buoyancy force, our system
can be modeled by the bond graph shown in Fig. 3 [15], [16].
The following forces are included: the wave excitation force
F
e
,
the inertia forces due to buoy mass and added mass at infinite
¨
¨
frequency
m
b
y
and
m
∞,b
y
. Furthermore, the hydrostatic stiff-
ness force
F
c
=
−C
b
y,
the damping force due to wave radiation
F
r
, and the load force from the power take-off machinery
F
l
is
accounted as well in the model.
As may be seen, the inertia force due to added mass at in-
finite frequency
m
∞,b
has been separated from the rest of the
radiation force
F
r
, in accordance with the usual time-domain
representation [5]. Taking the hydrostatic stiffness force to be
linear is equivalent to assuming a constant water-plane area
A
ω
for the buoy. The buoy is modeled as a sphere with radius
r
= 5[m].
The hydrostatic stiffness coefficient is then given by
C
b
=
ρgA
ω
. As long as the incoming waves are small com-
pared to the characteristic length of the buoy and the oscillation
amplitude is restricted to
|y| ≤
3
[m], this approximation keeps
force errors within
±6%
[17].
The bond graph has two energy-storing elements in integral
causality, which indicates that the model has two state variables:
the linear momentum
p
b
of the buoy, and its vertical position
y.
Also, we see that the force of the added mass is in differential
causality, which will provide us with an implicit equation for
the inertia force.
SCHOEN
et al.:
WAVE PREDICTION AND ROBUST CONTROL OF HEAVING WAVE ENERGY DEVICES FOR IRREGULAR WAVES
631
Fig. 2. (a) Influence coefficients of the predictive ARMA filter. (b) Influence coefficients of the predictive AR filter. (c) Influence coefficients of the proposed
predictive Kautz/AR filter. (d) Transformed filter weight distribution for the proposed predictive Kautz filter.
The hydrodynamic radiation force can be given as
F
r
=
t
−∞
κ(t
−
τ
)(p
b
/m
b
)dτ
, [5]. Here, the impulse response func-
tion
κ(t)
is a memory function linking the current radiation
force to past and current values for the buoy velocity
p
b
/m
b
.
The convolution integral
F
r
may be approximated by a finite-
order state-space model [18]
p
b
˙
ψ(t)
=
A
k
ψ(t)
+
B
k
m
b
Fig. 3. Bond graph for a heaving sphere showing the governing effects of
inertia, compliance, wave radiation, and excitation.
F
r
(t) =
C
k
ψ(t)
(16)
which will require to find the additional state vector
ψ.
Rear-
ranging (14), one arrives at
p
b
=
−
˙
m
b
(−C
b
y
−
F
l
−
F
r
+
F
e
).
m
b
+
m
∞,b
(17)
Deriving state equations from the bond graph yields
m
∞,b
p
b
−
C
b
y
−
F
l
−
F
r
+
F
e
˙
p
b
=
−
˙
m
b
p
b
.
y
=
˙
m
b
(14)
(15)
Defining the new state vector
x
= [
p
b
y ψ
T
],
one can now write
the entire system, including the radiation force approximation,
京公网安备 11010802033920号
Copyright © 2005-2024 EEWORLD.com.cn, Inc. All rights reserved
评论