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. Uncertainties in the physical description of the system as well as in the input from irregular waves provide challenges to the control algorithm. 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. INTRODUCTION 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 horizon 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 T Manuscript received May 29, 2010; revised September 1, 2010; accepted December 1, 2010. Date of publication February 4, 2011; date of current version 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 simulation results presented indicate reasonable success for long distances. One of the drawbacks of using the approach of Halliday 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 systems, 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 discussion of optimal control and phase control for WECs. Latching control attempts to achieve the same effect by halting the device’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 future 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 minima 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 prediction 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 environment 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. output vector. Defining Xss = ΦTAR YAR , the updated parameter estimate is θ̂AR (k + 1|k) = P (k + 1) × Xss (k + 1). Introducing the modal matrix Q composed of the eigenvectors of P = (ΦTAR ΦAR ) and diagonal eigenvalue matrix Λ of the correlation matrix P, the transformed filter weight vector results into II. WAVE PREDICTION MODELS ′ ′ XSS (k + 1) = QΛQT × Qθ̂AR (k + 1|k) = QΛθ̂AR (k + 1|k) ′ θ̂AR (k + 1|k) = QT θ̂AR (k + 1|k). (5) Transforming Xss : Pθ̂AR (k + 1|k) = Xss (k + 1) A. AR Predictive Model Consider a simple discrete time finite impulse response (FIR) AR model to characterize irregular waves ŷ(k|k − 1) = p  ′ ′ Xss (k + 1) = Λθ̂AR (k + 1|k). aj y (k − j) (1) j=1 where the output y ∈ Rn o x1 is the wave height with no = 1 and is measured up to time t = k − 1, k is the discrete time index, p is the model order, aj ∈ Rn 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 decomposition theorem, an AR model can fit a variety of data characteristics, 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 established (see [11]) yˆ(k + i|k − 1) = p  (i) aj 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: (i) (i−1) x′i (i−1) aj + a j + 1 . (3) 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 operation 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 property 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 = (ΦTAR ΦAR )−1 ΦTAR YAR (4) where θ̂AR is a vector containing the estimated parameter matrices, ΦAR is the information matrix, and YAR is the respective (6) x′ss Using as the elements in and λi as the corresponding eigenvalues found on the diagonal of Λ, each individual transformed filter weight contribution is given by (0) a′ i x′i , λi = for 1 ≤ i ≤ p. (7) B. ARMA Predictive Model Adding to the AR formulation a MA part enhances the flexibility of the AR model [12]. Consider a typical ARMA model in the output error formulation yˆ(k|k − 1) = p  aj y (k − j) + p  bm ε(k − m) (8) m =0 j=1 j=1 aj = a1 ′ Defining Xss = ΦTAR Xss , which premultiplied with QT yields where aj and bj are the ARMA filter matrices. The MA associated with the bj matrices is nonlinear due to the prediction error formulation. The model can be constructed through backsubstitution of the output into the original filter model formulation. An i-step ARMA predictor thus can be given as yˆ(k + i|k − 1) = p  (i) aj y (k − j) + p  b(i) ε(k − m) (9) m  m =1 j=1 where the predictive parameters for the past residual terms are defined as (i) (i−1) (0) aj aj = a 1 (i−1) + aj + 1 for j = 1, 2, 3, . . . and ak = 0, for k > p. The predictive parameters for the MA terms are defined as (i−1) b(i) m = a1 (i−1) b(0) m + bm + 1 for m = 1, 2, 3, . . . and bk = 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 propose to choose orthonormal filters Li (q) such that (1) becomes y (k + i|k − 1) = p  c̃i Li (q)y (k). (10) i= 1 To keep the selected OBF model linear in the parameters, Li (q) are chosen before the estimation of the filter coefficient matrices c̃i occurs. Li (q)’s are chosen with the purpose of incorporating prior knowledge of the system to be modeled. Since the waves exhibit an underdamped dynamical characteristic, we choose Kautz filters for the Li (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) With this, the filter Li (q) can now be computed using the following orthonormalization process:  (1 − c2 )(1 − d2 ) (12) L1 (q, c, d) = 2 q + c(d − 1)q − d  (1 − d2 )(q − c) L2 (q, c, d) = 2 (13) q + c(d − 1)q − d and the higher order filter parameters by using the recursive formula as follows: −dq 2 + c(d − 1)q + 1 L2i−1 (q, c, d) = L2(i−1)−1 (q, c, d) 2 q + c(d − 1)q − d The limits on the coefficients c and d are given as and − 1 < d < 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 measurements y(k), with a fixed sampling time. Normalize the wave height time series using y(k) = y(k)/ym ax , where ym ax = sup{y(k)}. ⎡ ⎢ ⎢ ΦKautz (q, c, d) = ⎢ ⎣ Influence of Kautz parameters on prediction. 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 c̃i using (i) Θ̂Kautz = {ΦTKautz ΦKautz }−1 ΦTKautz YKautz (M − i) 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 = [c̃1 , c̃2 , c̃3 , . . . , c̃p ]T YKautz = [y(p + i + 1), . . . , y(M + i)]T . −dq 2 + c(d − 1)q + 1 . L2i (q, c, d) = L2(i−1) (q, c, d) 2 q + c(d − 1)q − d −1 < c < 1 Fig. 1. 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 executing Step 2. For this, we use a set of simulations where the average correlation coefficient is computed between the predicted 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 L5 = 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 standard deviation (STD) for a varying c coefficient, corresponding L1 (q, c, d)y (p) L1 (q, c, d)y (p + 1) .. . L3 (q, c, d)y (p) L3 (q, c, d)y (p + 1) .. . ... ... .. . Lp (q, c, d)y (p) Lp (q, c, d)y (p + 1) .. . L1 (q, c, d)y (M − 1) L3 (q, c, d)y (M − 1) ... Lp (q, c, d)y (M − 1) ⎤ ⎥ ⎥ ⎥ ⎦ 630 IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011 TABLE I CORRELATION FOR DIFFERENT PREDICTIVE FILTERS UNDER VARYING NOISE ENVIRONMENT 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 relationship 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 prediction model, the damping, we are interested in observing the influence of this parameter on the predicted time series. Utilizing 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 perception 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 comparison 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. BUOY MODEL 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 Fe , the inertia forces due to buoy mass and added mass at infinite frequency mb ÿ and m∞,b ÿ. Furthermore, the hydrostatic stiffness force Fc = −Cb y, the damping force due to wave radiation Fr , and the load force from the power take-off machinery Fl is accounted as well in the model. As may be seen, the inertia force due to added mass at infinite frequency m∞,b has been separated from the rest of the radiation force Fr , 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 Cb = ρgAω . As long as the incoming waves are small compared 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 pb 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 Fr = t −∞ κ(t − τ )(pb /mb )dτ , [5]. Here, the impulse response function κ(t) is a memory function linking the current radiation force to past and current values for the buoy velocity pb /mb . The convolution integral Fr may be approximated by a finiteorder state-space model [18] ψ̇(t) = Ak ψ(t) + Bk Fig. 3. Bond graph for a heaving sphere showing the governing effects of inertia, compliance, wave radiation, and excitation. Deriving state equations from the bond graph yields m∞,b ṗb − Cb y − Fl − Fr + Fe ṗb = − mb pb . ẏ = mb (15) Fr (t) = Ck ψ(t) (16) which will require to find the additional state vector ψ. Rearranging (14), one arrives at ṗb = − (14) pb mb mb (−Cb y − Fl − Fr + Fe ). mb + m∞,b (17) Defining the new state vector x = [ pb y ψ T ], one can now write the entire system, including the radiation force approximation, 632 IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011 Fig. 4. Control system architecture. Dashed lines are associated with off-line operations. Fig. 5. FIS interaction with the WEC. as ⎡ −m b C b m b +m ∞, b 0 ⎢ ẋ = ⎣ m1 ⎢ +⎣ 0 Ak −m b m b +m ∞, b mb m b +m ∞, b 0 0 0 0 = Ãx + B̃ Fl Fe ⎤ ⎥ 0 ⎦x 0 b Bk ⎡ −Ck . ⎤ ⎥ Fl ⎦ Fe (18) In this form, Fl is the force from the PTO machinery on the buoy. This PTO force is the force used to control the system and is detailed in the following sections of this paper. For the studied heaving sphere, we have used data from [19] and [20] to compute the memory function κ(t). A reasonable approximation for the convolution term can be achieved with a state model of order six using Hankel singular value decomposition to reduce the model order [18]. This gives a total of eight states for the heaving buoy model, taking the wave excitation force and the machinery load force as inputs to the system. IV. FL CONTROL In this section, we propose to use a FL controller with the purpose to adapt the immediate PTO force Fl = −cẏ by controlling the damping coefficient c. The fuzzy controller includes a number of rules that allow for reducing the buoy motion, if the excursion becomes too large. The control output is computed based on the prediction of the future buoy motion, as a result of current measured wave heights. The horizon for this controller is of a rather short distance, i.e., 1 s. To accomplish a fuzzy inference system (FIS) for WEC control, the following system structure is proposed (see Fig. 4). In the figure, z denotes the wave height, y is the buoy motion, E is the energy, c∗ is the optimal damping for a given sea state, and the control parameter is the damping generated by the PTO. For the simulations, the wave energy device given by block WEC is defined by the model given in Section III [18]. The system identification SI in Fig. 4 provides a state-space model (resulting into the statespace matrices Ai , Bi , Ci , and Di ) to simulate the dynamics of the buoy for given predicted wave motions. The SI is performed off-line before the use of the controller. The prediction of the wave height is accomplished by a discrete time predictive AR formulation, which can be achieved by using (2). The proposed FIS controller implementation includes as the control variable the damping provided by the PTO. The Mamdani type FIS is defined for this control action by having three inputs and one defuzzyfied output, as depicted in Fig. 5: the three inputs are processed using a set of different input membership functions. For the Amplitude, the membership functions are categorized as “‘Small Difference,” “Some Difference,” “‘Large Difference,” and “‘Excessive Difference,” where Difference is defined as the discrepancy between |ŷ − z| and the maximum excursion allowed. The selected membership function for the amplitude input are sigmoid functions for the “Small” (left) and “excessive difference,” (right), while the remaining membership functions are based on the Bell membership function. A similar set of membership functions are created for the second input, the phase difference between the buoy motion and the wave motion. The phase is computed by establishing a time record of incoming waves and the buoy’s motion. The “Small Phase” and “Large Phase” are characterized with left sigmoid functions and right sigmoid functions, respectively. The “Good Phase” case is modeled with a Bell membership function. The last input, the velocity, representing the buoy’s vertical velocity, is used for the computation of the amount of damping. The single membership function for the velocity is given by a Polynomial-Z membership function. For the output, the damping is constructed by the use of four different membership SCHOEN et al.: WAVE PREDICTION AND ROBUST CONTROL OF HEAVING WAVE ENERGY DEVICES FOR IRREGULAR WAVES Fig. 6. Fig. 7. Buoy motion and wave motion under control. Fig. 8. Buoy motion with and without control using optimal parameters. 633 FS output surface for two inputs. functions, i.e., “Less Damping,” “Optimal Damping,” “More Damping,” and “Halting.” The corresponding membership functions are all Gaussian type membership function with the exception of the “Less Damping” case for which a sigmoid function is used, and the “Halting” case, where a PolynomialZ membership function is employed. The output damping is then computed as a change from the nominal value of damping. The nominal damping is equal to the optimal damping computed using a GA, which in turn is based on past input/output data, with the objective to absorb the maximum energy. The employed GA is based on a continuous number representation, with an elitism-based selection scheme. For more details on the used GA, see [21]. For the simulations, 100 generations are used with an initial population size of 96 chromosomes, and a steady-state population size of 48 chromosomes. The mutation rate is set to 4%. The input functions are linked to the controller output by a rule base of the form as follows: 1) if (Phase is Small_Phase), then (Damping is Less_ Damping); 2) if (Phase is Good_Phase), then (Damping is Optimal); 3) if (Phase is Large_Phase), then (Damping is More_ Damping); 4) if (Amplitude is Small_Difference), then (Damping is More_Damping); 5) if (Amplitude is Some_Difference), then (Damping is Optimal); 6) if (Amplitude is Large_Difference), then (Damping is More_Damping); 7) if (Amplitude is Excessive_Difference), then (Damping is Halting); 8) if (Velocity is zero), then (Damping is Halting). The overall resulting transfer output surface is given by Fig. 6. We recognize that the nominal amplitude difference and the phase is normalized to 0 (desired state). The defuzzyfication is done using the centroid of the output membership function. For the simulations, a Pierson–Moskowitz wave spectra of 46.3 km/h [25 knots, strong breeze] wind speed is utilized, while predicting the future wave height 10 steps (1 s) ahead. Fig. 7 depicts the simulated buoy motion under the influence of the FL controller, which is turned ON after 65 s (k = 650). At the same time of switching the controller ON, the optimal initial damping (and stiffness) is incorporated, which is used as the nominal damping from which the controller deviates. Fig. 8 provides a depiction of the buoy’s motion when optimal damping is used without control compared to with control action. The excursion limits are exceeded when no control is employed, while the proposed controller keeps this excursion within the desired limit of ±3[m]. The damping over time is depicted in Fig. 9. The short-term forecast and reaction by the FL controller is clearly manifested in the damping adaptation plot. The converted energy with and without the proposed FL controller is depicted in Fig. 10 [the energy is computed according to (26)]. Here in both cases, the optimal parameters are used. The energy graph indicates that the improvement is considerable. Although the corresponding buoy motions as depicted in Fig. 8 shows larger excursions for the uncontrolled case, the small change in the phase allows for a larger energy accumulation. It is also to 634 IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011 Fig. 11. Rules with corresponding membership functions and firing strength of output. Fig. 9. Damping computed using the FL controller. Fig. 12. Fig. 10. Energy generated with and without FL controller using optimal parameters. note that the membership functions of the inputs and the output are generated using simple rules and intuitive assessment of the consequences; hence, they are by no means optimal. The FL inference system detailed earlier is expanded and a second output is added to influence the stiffness of the PTO in addition to the damping. The PTO force then takes the form Fl = −cẏ − ky, where k is the stiffness coefficient to be controlled along with the damping coefficient c. In Fig. 4, the “Parameter” signal flow includes now also the stiffness in addition to the damping. The FL controller has two output signals, an addition of three rules, while maintaining the number of inputs. The three inputs are processed using a set of input membership functions. The output damping and stiffness are then computed as a change from the nominal values, which are the optimal coefficients computed using a GA based on past input/output data. The input functions are linked to the controller output by a rule base as depicted in Fig. 11. The computed optimal damping and Time history of buoy motion and wave motion. the optimal stiffness, using the aforementioned GA, are taken as the initial values when the controller is switched ON. After that, the controller adapts the new damping and stiffness for the immediate future only. Fig. 12 depicts the response of the buoy over time. The controller is turned ON at time k = 650. The immediate reaction of the device is to oscillate in resonance. This oscillation always stays within the excursion limits of ±3 m. Fig. 13 depicts one output of the controller, the stiffness (the damping plot is similar to Fig. 9). The converted energy is correspondingly increasing on a steeper curve than when no control is used (see Fig. 14). The optimized structural coefficients take on the following values: c∗ = (16.429e + 5) Ns/m and k ∗ = (1.0742e + 5) N/m, respectively. Also included in Fig. 14 is the comparison between the damping controller and the damping and stiffness controller. Fig. 14 indicates that a small improvement of the energy output is noticed due to the help of adapting the stiffness of the PTO during the simulations. The change is not significant, but indicates that an instantaneous update of the control parameters can have some influence. V. ROBUST CONTROL The proposed control architecture involves the identification of a generic WEC and the prediction of future wave height data. SCHOEN et al.: WAVE PREDICTION AND ROBUST CONTROL OF HEAVING WAVE ENERGY DEVICES FOR IRREGULAR WAVES 635 Elsner and He (in [23]) refined this as d (A, B) = min σn ([A − sI, B]) = min σ (s) s∈C Fig. 13. s∈C (20) where σn (G) is the nth singular value of the (n × (n × m)) matrix G. This implies that the problem of finding the distance to uncontrollability is the problem of minimizing σ(s) over the complex plane. It is shown in [23] that the function f (s) = vnH (s)[ un (s) 0 ]T is the partial derivative of σ(s) with respect to the components of s. The normalized left singular vector un (s) and the normalized right vector vn (s) are the nth column of U and V, respectively, which results from the singular value decomposition [A − sI, B] = UΣVH , and H denotes the complex (Hermitian) transpose. In [23] the function f(s) is used to find the minimum distance to the uncontrollable region by use of an iterative Newton method. The disadvantage of this method is that a good starting value for the iteration is needed and that no guaranty is given that the global optimum is found. We utilize an alternative approach by employing a simple GA. The search area can be defined ( [23]) by     A + AT A + AT ≤ x∗ ≤ λm ax Ireal = λm in 2 2     A − AT A − AT Iim ag = λm in ≤ y ∗ ≤ λm ax 2i 2i Stiffness (N/m) computed using the FL controller. (21) where s∗ = x∗ + iy ∗ defines the optimum point in the complex s-plane. To implement a robustness measure in the controller design, consider the closed-loop system given by (A + BF), where F ∈ Rm ∗ n is the feedback gain matrix. Assuming an uncertainty about the system reflected by ∆A and ∆B, then we need to maximize the distance d d(A + ∆A + (B + ∆B)F, S) Fig. 14. Comparison of Energy generated. Both processes include inaccuracies and errors that can lead to undesirable performance reduction of the WEC. Hence, we are interested in designing a controller that has some sort of robustness in terms of model inaccuracies. The objective of the controller proposed in this section is to limit the loss of robustness due to the closed-loop arrangement while maximizing the energy output of the WEC. The robustness is achieved with an optimization algorithm that maximizes the distance to the unstable region.Consider a continuous time, linear time-invariant system given by x˙ = Ax + Bu (19) where A ∈ Rn ∗n and B ∈ Rn ∗m . The system (A,B) is control¯ B̄]) = n ∀s ∈ C. lable, if rank([Ā − sI, Paige (in [22]) defines the distance to uncontrollability as the spectral norm distance of the pair (A,B) from the set of all uncontrollable pairs d(A, B) = min{||[E, H]|| : (A + E, B + H) uncontrollable} (22) where S ∈ Rn ∗n is the set of all unstable matrices. Assuming that A + BF + (∆A + ∆BF) has at least one eigenvalue with a positive real part, that is Re(λ) ≥ 0, we can define the robustness criteria as follows: the distance between the closed-loop system (A + BF) to the region of instability must be larger than the uncertainty in the closed-loop system. Expressed with aforementioned formulation (see [24]) d2 (A + BF, S) ≤ ||∆A + ∆BF||2 (23) where we made use of the spectral norm. Following the proposition of Kenney and Laub in [24], we can assert the following, For a closed-loop system (A + BF) with a system uncertainty of ∆A and ∆B, the spectral norm of the distance to the set of all unstable matrices S is given by d(A + BF, S) ≤ (1 + ||F||)d((A, B), S). Using (3), the above equation is equivalent to d(A + BF, S) = min σn ([Ac − λI]) Re( λ)≥0 where Ac = A + BF. Therefore, we propose to design the robust controller as follows: find F such that the minimum distance 636 Fig. 15. IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011 Controller architecture for optimal robust FL predictive control. Dashed lines are associated with off-line operations. to uncontrollability is maximized subjected to the constraint given in (21). Expressed in matrix form Max{Min{σn ([A + BF − sI, B])}} subjected to Min{σn ([A + BF − λI]) } ≤ (1 + ||F||) + Min{σn ([A − s∗ I, B])} and Ireal m i n ≤ x∗ ≤ Ireal m a x and Iim ag m i n ≤ y ∗ ≤ Iim ag m a x . The overall controller is achieved by combining the proposed robust controller with the proposed FL controller. In this fashion, we can address the short-term fluctuations of the waves by tuning the optimal damping force of the WEC device using the FL controller, and the long-term fluctuations are addressed by the GA, while inaccuracies in the model and the variations of conditions are in part compensated by the proposed robust controller. The resulting controller architecture/design is given in Fig. 15. The implementation of this controller design is done using a simple continuous GA (see [21]) that finds the appropriate minima for the defined cost function. The controller design for the robust control portion is based on a linear quadratic (LQ) controller. We have to pay special attention to the plant input, since only one input will be used to control the WEC, while both inputs are used to compute the control force. Considering the input to the WEC as u = [ z Fl ]T , where FL is the damping force, the state-space form of the system can be expressed as follows:     K1 K1 ẋ = A − B x and y = C − D x K2 K2 where we partitioned the feedback-gain matrix according to the input u. The controller will not have any influence on the wave force exerted onto the WEC; hence, K1 is set to zero and the Fig. 16. Accumulated absorbed energy using robust controller and no controller. TABLE II SIMULATION RESULTS FOR ACCUMULATED ABSORBED ENERGY FOR 200 SECONDS, WITH AND WITHOUT HYBRID ROBUST FL CONTROL closed-loop system representation can be given as ẋ = [A − B1 K2 ]x y = [C − D2 K2 ]x with B = [ B1 B2 ] and D = [ D1 D1 ]. Provided that u2 is sufficient to reach all system states, the updated optimization becomes: Find K2 such that SCHOEN et al.: WAVE PREDICTION AND ROBUST CONTROL OF HEAVING WAVE ENERGY DEVICES FOR IRREGULAR WAVES 637 where the output error model is given as y(k) = Fig. 17. Buoy motion with and without robust FL controller. J1 = Max{Min{σn ([A + B1 K2 − sI, B1 ])}} (24) subjected to Min{σn ([A + B1 K2 − λI])} ≤ (1 + ||K2 ||) + Min{σn ([A − s∗ I, B1 ])} and Ireal m i n ≤ x∗ ≤ Ireal m a x and Iim ag m i n ≤ y ∗ ≤ Iim ag m a x . In the following, the introduced robust controller is implemented in discrete time. We use the traditional discrete time LQR controller design for K2 , which is given as follows:   JLQ R = xT (k)Qx(k) + uT (k)Ru(k) (25) where Q and R are used as the variables in the GA in order to optimize the discrete time version of (24). The energy production is given by the following cost function:  ẏ T (k)cPTO (k)ẏ(k). (26) JE = The overall optimization can be expressed by: maximize the minimum singular value of the closed-loop system (24) while maximizing the energy output of the WEC (26) using as little control energy as possible (25). The identification of the model is accomplished once at the start of the operation of the WEC and may be repeated after some time to update the system model. The identification is based on the actuation force and the wave as inputs, where the actuation force is varied randomly with increasing amplitude according to F (k) = {10 × 1.025k } × η, where η = N (1, 0). The discrete time model can be given as an output error model with the following polynomials: B(q) u(k) + e(k) H(q) (27) and e(k) is the residual vector. The proposed controller is implemented in two phases to investigate each component, i.e., the robustness issue and the short-term interaction of the FL controller. Fig. 16 depicts one simulation run for the robust controller portion in terms of energy absorbed. In the simulation, the robustness is tested using a perturbation of the system model characteristics, i.e., one of the parameters of the buoy model is altered. In example, the WEC model assumes that the hydrostatic stiffness force is linear, which implies a constant waterplane area for the buoy. Since this is only an approximation, we use this parameter input by varying it from its nominal value. The controller is switched ON at t = 20 s, after completing the identification task to extract the model. For the simulation duration and selected sea state, the controller enables the system to convert more energy than the uncontrolled case (without altering the parameters on a continuous base). Also to note is that the control action is powered from the stored energy; hence, control forces aligned with the buoy velocity reduce accumulated absorbed energy by the amount of the control energy, while opposed control forces extract energy. Perturbing the system as described earlier, the simulations are repeated and the influence of the controller is compared to the unperturbed case. The resulting energy production loss due to the change (6%) in the selected model parameter is 19% for the controlled case, and 77% for the uncontrolled case. Combining the robust controller and the FL controller, as outlined in Fig. 15, the uncertainty influence can be quantified through simulations. The optimal stiffness and damping values computed with a GA for a given sea state are used as the nominal values, while perturbing the resulting hydrostatic stiffness force as the uncertainty in the WEC model and regulating the damping through the proposed controllers. Table II summarizes the results for a varying perturbation of the hydrostatic stiffness force. Table II reveals the influence of the proposed robust FL controller to changes in absorbed energy by the WEC with respect to uncertainty in the hydrostatic stiffness force. The controller has the effect to minimize the affect such an uncertainty may have in the operation and ultimately the energy production of a WEC. Fig. 17 depicts the comparison of the buoy motion for the case of using the combined (hybrid) robust FL controller (turned ON at k = 650) and no control for nonoptimized parameters (stiffness and damping). The small changes indicated by differences in the amplitude and the phase between the controlled and uncontrolled case cause substantial differences in the energy conversion as indicated by Fig. 16 extrapolated for long durations of operations. B1 (q) = 0.05616q −1 − 0.04712q −2 B2 (q) = −4.523 × 10−8 q −1 + 4.518 × 10−8 q −2 H1 (q) = 1 − 1.988q −1 + 0.9893q −2 H2 (q) = 1 − 2q −1 + 1.002q −2 VI. CONCLUSION The two proposed prediction algorithms (Kautz/AR and the Kautz) indicate to perform in a competitive fashion when compared to existing prediction models. The reduced model or- 638 IEEE TRANSACTIONS ON ENERGY CONVERSION, VOL. 26, NO. 2, JUNE 2011 der of the proposed prediction filters will help lower the rather high computational load imposed by conventional filters. For the operation of WECs in irregular waves, two controllers are proposed, each having a different objective. The FL controller addresses the short-term control, while the robust controller attempts to minimize the effects of the fluctuations or errors of the WEC model parameters on the energy absorption. The control architecture proposed requires the prediction of wave height data as well as the optimization of the structural parameters of the vibrational system, which is accomplished with a simple GA. Identification of the model is achieved off-line and may be repeated after some time when wave conditions have changed. Simulation results for each controller are presented separately before the combined robust FL controller is tested. The proposed combined controller has the effect to limit the excursion of the buoy motion to the predesigned limits, while maximizing the energy output. In addition, uncertainties in the model have a more limited effect on the energy production when the proposed controller is used. REFERENCES [1] J. L. Forsberg, “On-line identification and prediction of waves,” in Hydrodynamics of Ocean-Wave Energy Utilization, D. V Evans and A. F. de O Falcao, Eds. Berlin, Germany: Springer-Verlag, 1986, pp. 185–193. [2] J. R. Halliday, D. G. Dorell, and A. Wood, “A Fourier approach to short term wave prediction,” in Proc. 16th Int. Offshore Polar Eng. Conf., San Francisco, CA, 2006, pp. 444–451. [3] R. H. Jones, “Identification and autoregressive spectrum estimation,” IEEE Trans. Autom. Control, vol. AC-19, no. 6, pp. 894–897, Dec. 1974. [4] L. Budal and J. Falnes, “A resonant point absorber of ocean-wave power,” Nature, vol. 256, pp. 478–479, 1975. [5] J. Falnes, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave Energy Extraction. Cambridge, U.K.: Cambridge Univ. Press, 2002. [6] U. A. Korde, “Efficient primary energy conversion in irregular waves,” Ocean Eng., vol. 26, pp. 625–651, 1999. [7] U. A. Korde, “Latching control of deep water wave energy devices using an active reference,” Ocean Eng., vol. 29, pp. 1343–1355, 2002. [8] A. Babarit and A. H. Clement, “Optimal latching control of a wave energy device in regular and irregular waves,” Appl. Ocean Res., vol. 28, pp. 77– 91, 2006. [9] A. Babarit, G. Cuclos, and A. H. Clement, “Benefit of latching control for a heaving wave energy device in random sea,” in Proc. 2003 Int. Offshore Polar Eng. Conf., Honolulu, Hawaii, 2003, pp. 341–348. [10] S. M. Key and S. L. Marple, “Spectral analysis—A modern perspective,” in Proc. IEEE, Nov., 1981, vol. 69, no. 11, pp. 1380–1419. [11] U. A. Korde, M. P. Schoen, F. Lin, “Time domain control of a singlemode wave energy device,” Presented at the Int. Soc. Offshore Polar Eng., Stavanger, Norway, 2001. [12] O. Nelles, Nonlinear System Identification from Classical Approaches to Neural Networks and Fuzzy Models. Berlin, Germany: Springer, 2001. [13] P. J. Brockwell and R. A. Davis, Introduction to Time Series Analysis and Forecasting. New York: Springer, 1997. [14] W. J. Pierson and L. A. Moskowitz, “Proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii,” J. Geophys. Res., vol. 69, pp. 5181–5190, 1964. [15] D. C. Karnopp, D. L. Margolis, and R. C. Rosenberg, System Dynamics: Modeling and Simulation of Mechatronic Systems, 4th ed. New York: Wiley, 2006. [16] H. Engja and J. Hals, “Modeling and simulation of sea wave power conversion systems,” presented at the Eur. Wave Tidal Energy Conf., Porto, Portugal, 2007. [17] J. Hals, T. Bjarte-Larsson, and J. Falnes, “Optimum reactive control and control by latching of a wave-absorbing semisubmerged heaving sphere,” in Proc. Int. Conf. Offshore Mech. Arctic Eng., 2002, pp. 415–423. [18] R. Taghipour, T. Perez, and T. Moan, “Hybrid frequency—Time domain models for dynamic response analysis of marine structures,” Ocean Eng., vol. 35, pp. 685–705, 2008. [19] T. Havelock, “Waves due to a floating sphere making periodic heaving oscillations,” in Proc. Royal Soc., vol. 231 A, pp. 1–7, 1955. [20] A. Hulme, “The wave forces acting on a floating hemisphere undergoing forced periodic oscillations,” J. Fluid Mech., vol. 121, pp. 443–463, 1982. [21] R. L. Haupt and S. E. Haupt, Practical Genetic Algorithms. New York: Wiley, 1998. [22] C. C. Paige, “Properties of numerical algorithms relating to controllability,” IEEE Trans. Automat. Control, vol. AC-26, no. 1, pp. 130–138, Feb. 1981. [23] L. Elsner and C. He, “An algorithm for computing the distance to uncontrollability,” Syst. Control Lett., vol. 17, pp. 453–464, 1991. [24] C. Kenney and A. J. Laub, “Controllability and stability radii for companion form systems,” Math. Contr. Signals Syst., vol. 1, pp. 239–256, 1988. Marco P. Schoen (SM’10) was born in 1965. He received the B.S. degree in mechanical engineering from the Swiss College of Engineering, in 1989, the M.E. degree in mechanical engineering from Widener University, Chester, PA, in 1993, and the Ph.D. degree in Engineering Mechanics from Old Dominion University, Norfolk, VA, in 1997. From 1997 to 1998, he was a Faculty Member at Lake Superior State University, and from 1998 to 2001, he was a Faculty of the Mechanical Engineering program at Indiana Institute of Technology. Since 2001, he has been with Idaho State University, Pocatello, where he is currently a Professor and Chair in the Department of Mechanical Engineering, and an Associate Director of the Measurement and Control Engineering Research Center (MCERC). He has been an Associate Editor of the Journal of Dynamic Systems, Measurement and Control. His research interests include controls and vibration of biomedical and aerospace systems as well as energy-related problems. Dr. Schoen has been the Past Chair of the Model Identification and Intelligent Systems (MIIS) Technical Committee for the American Society of Mechanical Engineers (ASME). Jørgen Hals was born in 1974. He received the M.Sc. degree in physics from the Norwegian University of Science and Technology (NTNU), Trondheim, Norway, in 1999, and the Ph.D. degree from Norwegian University of Science and Technology (NTNU), in 2010. His Ph.D. study was focused on renewable energy (wave energy conversion). He was a Research Fellow at NTNU and a Researcher at SINTEF Energy Research, where his work was focused on a wide range of energy technologies. He is currently a Postdoctoral Fellow at the Center for Ships and Ocean Structures (CeSOS), NTNU. His research interests include mathematical modeling and optimal control of wave energy converters. Torgeir Moan was born in 1944. He received the M.Sc. and Ph.D. degrees, in 1968 and 1976, respectively, from the Department of Civil Engineering, Norwegian University of Science and Technology (NTNU), Trondheim, Norway. Since 1977, he has been a Professor of marine structures at NTNU, where he has been the Director of the Research Center of Excellence (CeSOS), since 2002. He was a Visiting Professor at Massachusetts Institute of Technology for one year, at University of California, Berkeley for two years, and a part-time Adjunct (Keppel) Professor at the National University of Singapore during 2003–2008. He has authored more than 400 scientific papers, delivered more than 20 keynote, plenary lectures in international conferences and award lectures, as well as educated 50 doctoral candidates. He is an Editor of the Journal of Marine Structures and is on the Editorial Board of eight other journals. His research interests include stochastic dynamic modeling and reliability and analysis (relating to all kinds of marine structures), and facilities for offshore wind and wave energy. Dr. Moan has been an Elected Member of all the national academies of science and letters in Norway and a Foreign Fellow of the Royal Academy of Engineering, U.K., as well as a Fellow of several international professional societies. He has received numerous research prizes.