WO2013030164A2

WO2013030164A2 - Wave energy converters

Info

Publication number: WO2013030164A2
Application number: PCT/EP2012/066624
Authority: WO
Inventors: Giorgio BACELLI; John Ringwood; Jean-Christophe GILLOTEAUX
Original assignee: National University Of Ireland
Priority date: 2011-08-26
Filing date: 2012-08-27
Publication date: 2013-03-07
Also published as: WO2013030164A3

Abstract

By approximating the motion of a wave energy converter (WEC) device and of the excitation forces acting on the device, each as a finite, truncated combination of basis functions, the equations of motion describing the device can be reformulated as a linear system and the problem of the constrained energy maximization can be reformulated as a non linear program expressed in terms of a cost function and a constraint. This can be solved in real time to arrive at a force profile for a PTO which, when implemented by the PTO, maximises the energy available from the PTO while staying within the physical constraints imposed by the device.

Description

Wave Energy Converters

Technical Field

This invention relates to wave energy converters and in particular to methods, apparatuses and computer programs for controlling wave energy converters (WECs).

Background Art

Wave Energy Converters (WECs) are devices designed to extract energy from water waves.

Several technologies are under development, such as oscillating body, oscillating water column and overtopping device; the main distinguishing characteristic between them is the energy conversion process. The oscillating body, the type of device considered in the present case, recovers energy form the motion of floating bodies subject to the action of waves, by means of a mechanical device named Power Take Off (PTO). The PTO applies a force on the oscillating body, and the mechanical work performed by this force is the energy absorbed by the PTO, part of which will be converted to electricity and delivered to the grid. The oscillating body category of WECs can in turn be divided into single-body devices and self-reacting devices. For single-body devices, the PTO applies a force between the floating body and a fixed reference, such as the seabed. Self-reacting WECs are devices composed of multiple bodies, and energy is recovered by means of the force applied by the PTO between the bodies. The energy exchanged in the interaction between water waves and a floating body depends, among other factors, on the motion of the same body and on the motion of other floating bodies located in its proximity. Part of the energy carried by waves is transferred to the body by means of the excitation force, which is the force that waves exerts on the body. On the other hand, the floating body itself returns energy to water by radiating waves generated by its own motion. The radiated waves, in turn, affect the motion of any other floating body located in proximity and vice versa.

The force exerted by the PTO affects the motion of the device and, as a consequence, also influences the amount of energy exchanged between the water and the WEC. The two most common types of PTOs considered in the design of WECs are direct drive linear electrical generators and hydraulic cylinders connected to conventional rotational generators through hydraulic circuits. Both types of PTOs have limited ranges of operation due to the intrinsic physical limitations of the components used to build them. For example, hydraulic cylinders have a finite stroke, which is the difference between the total length of the cylinder when the rod is fully extended and the length of the cylinder when the rod is fully contracted. A self-reacting WEC composed of two bodies and equipped with an hydraulic cylinder as PTO is characterized by a limited displacement between the two bodies, that depends on the stroke of the cylinder. Besides, the maximum value for the oil pressure inside the cylinder specified by the constructor for safety operation limits the maximum force that the PTO can exert.

The benefits provided by the implementation of a control system to modify the dynamical behaviour of the WEC are widely studied in literature. The most significant improvement is the amount of energy produced, for example using the well know optimal control theory. Although it provides the mathematical means for the maximization of the extracted energy, it is of impractical implementation in real devices due to the requirement of future knowledge of the incoming waves. Besides, it does not provide tool to deal with motion and force restrictions introduced by the PTO. The same considerations apply for other control techniques studied in literature, except for Model Predictive Control (MPC), that has been successfully applied in simulation to a single-body device subject to motion and force restrictions. A significant drawback of MPC is its computational requirements, due to the large dimension of the non linear program that needs to be solved.

Disclosure of the Invention

The invention provides a method of controlling a wave energy converter as claimed in claim 1.

By approximating the motion of the device and of the excitation forces acting on the device, each as a finite, truncated combination of basis functions, it has been found that the equations of motion describing the device can be reformulated as a linear system and the problem of the constrained energy maximization can be reformulated as a non linear program expressed in terms of a cost function and a constraint. This can be solved in real time to arrive at a force profile for a PTO which, when implemented by the PTO, maximises the energy available from the PTO while staying within the physical constraint imposed. Preferred features of the invention can be found both in the dependent claims and in the particular description of the generalised solution which follows below, as well as in a more specific example based on the use of Fourier series as basis functions. It is to be understood that aspects of the invention are to be found both in the generalised equations and examples, and in the more particular Fourier- specific equations and examples.

A particular advantage of using a finite truncated combination of basis functions is that suitably selected basis functions allow modelling of the motion of the bodies in terms of a relatively restricted range of frequencies (for example, by selecting only the first N terms of a Fourier series). This restriction enables real time solutions to be found for the optimisation problem, without any significant inaccuracy in the solution due to the fact that the model which is thereby employed discards higher frequencies which will not in fact occur in a real-life, damped WEC with a PTO operating in waves found in nature. In other words, above a certain frequency, oscillation will not be found in the device, or at least not at any level which will significantly impact the accuracy of the real-time control imposed on the WEC.

There is also provided a computer program, a computing system, and a wave energy converter as claimed.

The computer program may be recorded on any suitable physical data carrier, or it may be transmitted as a signal. It may be hard-wired in a customised electronic circuit, or programmed into any suitable processor to provide a programmed apparatus which implements the method of the invention.

In a preferred embodiment, the control algorithm computes an approximation of the optimal motion of the device and of the optimal profile for the PTO force, based on the estimation and prediction of the excitation forces obtained from sensors placed on board of the device.

For optimal operation, it is intended the motion and the force profile that allows the maximum amount of energy to be converted without violating the force and motion restrictions. The sensors used for the estimation of the excitation force are preferably accelerometers and a sensor that provides the displacement of the PTO.

Brief Description of the Drawings

Fig. 1 is a block diagram of a system which employs the method of the invention;

Fig. 2 is a schematic of a self-absorbing point absorber system; Fig. 3 is a graph of relative position of a waverider buoy in a simulated system over time showing the effect of activating a control system and method of the invention in the period from t=100s to t=380s; and

Fig. 4 is a flowchart of operation of a method according to the invention.

Detailed Description of Preferred Embodiments

Fig. 1 is a block diagram of a system which performs the method of the invention. A wave energy converter (WEC) 10 is provided with a power take-off (PTO) 12 in known manner. A dedicated or general purpose computing system 16, which may be either local to or remote from the WEC 10, is programmed to implement a control algorithm. The computing system 16 receives inputs

14 representing the position (x) and acceleration (3 ) of the bodies making up the system. In the illustrated system the WEC is a two body system and thus the position can be the relative position of one body relative to the other, possibly as determined from the state of the PTO. The "Fe Pred. and Est." block 18 is a state estimator that reads the accelerations and the relative positions of the bodies, and it provides an estimate of the excitation forces. The predictor is implemented using an Auto Regressive model, as described further below in the detailed mathematical treatment. The "Reference Generator" block 20 computes the optimal motion and the optimal PTO force profile, based on the estimation and prediction of the excitation forces calculated in the block "Fe Pred. and Est." 18.

The "Reference Generator" provides an approximation of the motion and the PTO force that maximizes the energy absorbed by the device satisfying the restrictions that characterize the WEC. The estimated excitation force is first approximated by a linear combination of simple functions called basis functions, such as piecewise polynomial or trigonometric functions. The

approximated excitation force is then substituted into the mathematical model of the device together with the approximation of the motion and of the PTO force, using either the same or different basis functions. These substitutions provide the formula for the cost function that describes the amount of absorbed energy, to be expressed as function of the parameters describing the PTO force. The substitutions also provide the constraints to be expressed as function of the parameters describing the PTO. Optimal motion is then computed solving an optimization problem, with a dimension that depends on the number of basis functions used to describe the PTO force. Details of these calculations, both for the general case of basis functions in general, and for the specific case of a Fourier implementation, are provided in the detailed mathematical treatment below.

The "Low Level Controller" block 22 generates a control signal that makes the device follow the optimal trajectories computed by the "Reference Generator" 20. The control signal is applied to the PTO in order to exert the force required to allow the optimal motion.

The "Low Level Controller" has been designed as an optimal tracking problem with a reference also on the input: it tracks both position and velocity of the body (outputs), and also the PTO force, that is the control input. The tracking of the control input provides an increase of the absorbed energy.

The system provides the following benefits:

^■ Has the potential to double energy converted, compared to constant load PTO.

^■ Can naturally take into account the physical limitations of the device and PTO, such as motion restrictions and maximum PTO force.

■ The algorithm can be implemented in real time using information provided by easily available on-board sensors.

^■ A number of design parameters can be chosen by the user: type and number of the basis functions used for the computation of the optimal motion and PTO force; updating rate of the optimization; sensitivity of the optimization; type of predictor; prediction horizon.

■ The algorithm can take into account some non-linearities of the device.

^■ The algorithm can be implemented for the control of a single body WECs, single degree of freedom device or multi degree of freedom devices.

^■ The algorithm can be implemented for the control of multi-body WECs

^■ The algorithm can be implemented for the centralized control of an array of WECs. Detailed mathematical treatment explaining the justification for the invention, both at a generalised level and in a specific Fourier-based implementation, together with results obtained in a simulation. 1 Introduction

A Wave Energy Converter (WEC) is a device that produces electricity by converting the energy carried by water waves, usually ocean waves. The WEC considered in this work is a self-reacting point absorber; a floating body is said to be a point absorber when its horizontal dimensions are small compared to the length of the incident wave. A self- reacting point absorber is a WEC composed of several point absorbers that converts energy from the relative motion between the bodies. The energy is recovered by a Power Take Off (PTO) unit, which is a mechanical device capable of exerting a force between the bodies of the WEC; the amount of energy flowing through the device can be controlled by acting on the force exerted by the PTO. The maximization of the energy produced by a two-body point absorber has been addressed in Korde [2003] where the problem is approached in the frequency domain. Mixed time- frequency domain analysis has been carried out in Hal et al. [2007] and in Falcao et al [2009] where two-body WECs equipped with hydraulic PTO are considered. In practice most wave energy converters are characterized by physical limitations, such as the maximum relative motion between the bodies; in Falnes [1999] the author studied the effect of motion restriction on the energy produced by a two-body point absorber. Optimal control problems for WEC subject to amplitude restriction have also been addressed for single body devices as in Korde [2001 ], Eidsmoen [1996a] and Eidsmoen [1996b]. In Hals [2010] the author applied model predictive control for the maximization of the energy produced and considered also a constraint on the maximum force that can be exerted by the PTO.

The focus of this paper is to present a method for the solution of the energy maximization problem for a self-reacting point absorber subject to amplitude restriction. The objective is to provide an approximated solution of the constrained optimal control problem that requires a computational effort compatible with real time implementation. The problem is approached by approximating the force exerted by the PTO and the motion of the device by a linear combination of basis functions. The optimal control problem is then transformed into a constrained finite dimensional optimization problem.

The general formulation of the method is presented in section 2, while in section 3 an example is presented where the PTO force and the motion are approximated by means of Fourier series. Section 4 describes an optimization algorithm used to solve the constrained optimization problem; the practical implementation of the control system and simulation results are provided in section 5 and 6, respectively.

2 General formulation

The point absorber considered is composed of two concentric and axisymmetric bodies, both oscillating in heave (Fig. 2). Body A is a torus with rectangular section and mass m^A while body B is a cylinder with the axis directed along the vertical direction and with mass m^B . The vertical velocities of body A and body B are denoted v^A (t) and v^B{t) , respectively. Energy is recovered from the relative motion between the two bodies by means of the PTO, which is capable of applying a force f ,₀{t) between body A and body B . The energy absorbed by the PTO in a time interval of length T , J{T) , is equal to the work performed by f_plo(t) in the same interval, that is: T) = [ ,„( ( "( - ^Β( )<Λ,

thus, the energy absorbed by the device can be controlled via . The force exerted by the PTO also affects the motion of the bodies, which is described, for small oscillations, by the linear model (Falnes [1999]):

L^A (t) = m^Av^A (t) + B^Av^A (.t) + S^Au^A (t) - ff (0 - ff (/) + f_plo (0 = 0

= m^Bv^B(t) + B^Bv^B(t) + S^Bu^B(t) - / (/) - /_r ^fl( - f_plo(t) = 0 where u ^A(t) and u^B(t) are the vertical positions of the two bodies. The hydrostatic buoyancy is described by S^A and S^B while B^A and B^B are terms describing the linear viscous loss. The excitation forces on body A and on body B are denoted by f_e ^A and f , respectively. The excitation force is the force acting on a body due to the incident wave when the body is held fixed. The radiation forces f_r ^A and are forces acting on a body due to the radiated wave resulting from its own oscillation or by the oscillation of a second body located in proximity. The radiation forces may be expressed as functions of v^A and v" , evidencing the interaction between the bodies and the coupling of the equations in (2), as: m AA . A _{t) _ _kAA _(/)„ _VA _w _ _m B - B _(/) _ ^B ₍₍₎ * _yB _{()

m B . B _(/) _ _kBB _(t) , _vB ₍ _ _m BA . A _(/) _ £ BA _(/) , _yA _(/) where the symbol * denotes the convolution operator and the parameters w' and k^iJ (t) , with i,j e {Α, Β} , are the asymptotic values of the added mass and the radiation impulse responses, respectively.

The objective of the work presented in this paper can now be formally stated as: given the WEC described by the model in (2), find the optimal profile of the PTO force (f_plo ) in a given time interval of length T that maximizes the absorbed energy J(T) as defined in (1), subject to:

where the infinity norm is defined in (A. l). The constraint described by (3) is a restriction on the maximum distance between the two bodies and it is generally due to the PTO. For example, if the PTO is an electric linear generator or an hydraulic ram, which are the most common types of PTOs considered for this kind of WECs, the constraint defined in (3) may refer to the maximum excursion allowed by the linear generator or to the stroke length of the hydraulic piston.

The control of a two-body point absorber subject to motion restriction has already been addressed in Falnes [1999], where the problem for the two-body device is reformulated as the control of an equivalent WEC composed of only one body. The author also analyzed the effect of motion constraints on the average absorbed power, extending the work in Evans [1981 ] to a two-body device; the analysis is carried out in the frequency domain and the motion restriction is defined as Γ ||ΔΚ|| < ΔΚ„² _∞

with ||ΔΓ|| = AV'AV and A V^* denoting the conjugate transpose of AV . Equation (4) states a condition on the relative velocity AV between the bodies, but it can be easily converted to a condition on the relative position AU considering that V = iu> U . If the relative motion of the device is a sinusoidal oscillation, then the constraint ||Δί/|| < AO_m ² _ax is equivalent to (3). However, if the motion is still periodic with period T but not sinusoidal, the satisfaction of the frequency domain condition ||Δί/|| < AU²^ does not imply the satisfaction of (3); in fact, applying Parseval's theorem and (A.3) immediately shows that \\AU\\≤ ||w^( - «^s( |² · Thus, there may exist some values of t for which I u^A {t) - u^B {t) |> AU_max even if (4) is satisfied. Therefore, the formulation in (4) for the amplitude restriction is not adequate because (3) describes a physical limit of the device that cannot be violated.

The solution presented in this paper is based on the discretization, in the time domain, of the PTO force and of the motion of the device in order to transform the problem into a Non Linear Program (NLP). The approach is similar to the direct simultaneous method used for the solution of optimal control problems (Cuthrell and Biegler [1987]), where both the control variables and the state variables are discretized.

The PTO force is assumed to be such that f_plo (t) e L² ([0,T]) , where L² ([ , T ) is the Hilbert space of square integrable functions in the interval [Ο, Γ] ; also

v^A {t), v^B{t) e L² ([Q,T ) because they are velocities of physical bodies. The PTO force and the velocities are then approximated as a linear combination of basis functions in a finite dimensional subspace of the space L² ([0, T]) : ν^,( * "( =∑^Φ_>(

/,.^{(/) «} * ( =∑ Λ' ( where {φ, (/),... ,<lv(t)} is a basis for the finite dimensional subspace S^v c L²([0,T]) and is a basis for the finite dimensional subspace S^p a L² ([0,T]) . For any

given set of coefficients describing the PTO force {ρ. , . , . , ρ „} , the components of the

N

velocities are calculated by solving the system: ί ⁱ <^\ι-⁴( ,,Φr,,>^/ = ο V/ = 1, ... N

[<^( ,φ,> = 0

where (·,·) denotes the inner product defined as

Developing the calculations (Appendix B), the system (8) results in the linear system

with X^A

, and G is an N -by- N matrix, T is an N -by- N^p matrix and K is a vector of dimension N . The existence of a solution for the linear system (10) depends on the choice of the basis functions φ .(f) and on the hydrodynamic parameters characterizing the WEC, which determine the singularity of the matrix G . If a solution exists then, by an appropriate partitioning of the matrix G^~lT as [T^AT^BY and of the matrix G^'lK as [E^A E^Bf , the energy converted by the device, described by (1), can be written as a quadratic equation in the vector P :

J(P) = P^T(X^A - X^B) = P^TAT^TAP + AEAP

with = T - T^B , AE = E^A - E^B and the elements of the matrix A defined as

The approximated absolute vertical positions of the bodies are obtained by integrating (5) and (6):

and with u* = u^A (0) and u = u^B ( ) . With the substitution of u^A (t) and u^B {t) into (3), the constraint can be expressed as a function of the components of each velocity X^A and

X^B), which are related to the PTO force components through the affine transformation

(10). The result is that the constraint (3) becomes a function of the PTO force components P :

+ Φί ίΔΓ + ΔΕΐ≤AU_max, (12) where A ₀ = u^A - u and Φ(/) = [φ, (/),... ,φ_Ν (/)] .

The problem to be solved is now a finite dimensional constrained optimization problem in the variable P with the quadratic cost function in (1 1 ) and the inequality constraint in

(12).

3 Discretization using Fourier Series

Given the oscillatory nature of the problem, a truncated Fourier series is an intuitive

choice for the discretization and approximation of the problem. Furthermore, choosing

w₀ = 2π/Τ , the set of functions {sin (co_or) , cos(ro₀7') , ..., sin(Nro_o ) , cos(Nco_or)} form an orthogonal basis for the space with the inner product (9), thus giving the linear

system (10) and the cost function (1 1) favorable properties. The constant term of the basis is not considered because it is assumed that all the functions have zero mean; in practice, it is assumed that the reference frames of the bodies are chosen such that the origins

oscillate around their mean position with respect to the inertial reference frame. Using a zero mean truncated Fourier series with N frequency components for both the velocities and the PTO force, the dimension of each of the spaces S^v and S^p is 2N , and the resulting approximating functions (5)-(7) become:

N

v^A (t) = p_n ^A cos( >y) + sin(M(o₀ (13) n=l v^B ( =∑α_η ^Β cos(Hay) + b sin(mo₀ (14) n=l

N

f_plo(i) =∑α cos(n(D₀t) + b^B sin(«ro₀ . (15) For the practical implementation of the algorithm, described in section 5, it is convenient to approximate the excitation forces by a truncated Fourier series containing N frequency components: fe^A ( * fe^A ( = e cos(na>₀0 + e_B' sin(«m₀

n=l f (t) * f i = Y _n ^c cos(iK0₀i) + e 5ΐη(«ω₀ ·

The mean value of the excitation forces can be considered as zero with no loss of generality. In fact, the excitation force is calculated as the convolution of the wave elevation with the excitation force kernel (Falnes [1999]); the wave elevation can be transformed into a zero mean function by changing the origin of the reference frame, resulting in a zero mean excitation force.

3.1 Equation of motion

Substituting (13)-(17) into (2), the linear system corresponding to the equation of motion (10) is:

where I_2N is the identity matrix of size 2N .

X^A = [ ,b , a^A, b₂ , ... ,a^A ,b^A]^T,

rA iT

'-' L J :

p

S_\ , e₂ ' ' · · · ' ' 1 '

0 0 0

0 0 D^ij M

0 0 with i,j = {A, B} , and

D^ = R^AA (n( ₀) + B^A ,

M ^A = mo₀(m^A + (ηω₀))- S^A/(na)₀),

D^BB = R^BB (^₀) + B^B,

M?^B = «co₀(w^fl + ηι^ΒΒ(ηω₀))-

D^AB = R^AB(n *₀\ M^AB = ηω₀(ηι^ΑΒ(_ηωο)}

D^BA = R^BA (ηω₀), „^ω = «ω₀ (ι»¹" (ηω₀))

The radiation resistances i?'^v (»ro₀) and the added masses m'^J (no3₀) are related to the

impulse responses k'^j(t) and to the asymptotic added masses m' by means of the

equations in Ogi lvie [ 1964] .

The matrix G_v is block diagonal and each block is a 2-by-2 normal matrix of the form a b

(19) -b a

This particular structure is due to the orthogonality of the Fourier series and it allows the study of the existence of the solution of the linear system (18) by studying the singularity of each of the N 4-by-4 matrices

Each matrix G_n corresponds to a frequency mo₀ ; thus, should the system in (18) be

singular, a possible solution might be to perform a different frequency discretization by selecting a different fundamental frequency ω₀ .

3.2 Cost function

If the solution of (18) exists, the cost function (1 1 ) describing the amount of energy

absorbed by the PTO is

J(P) = -P^THP + P^T (Q^AE^A (21) where

H = S;_BB +

-1

Q ^ = G G^A"S G^~ A* A + ' S " GrAA -

S _M and S _BB are, respectively, the Schur complements of G and G , and they are defined as

It can be shown that the matrix H is positive definite because it is block diagonal with 2- by-2 blocks of the form (19); therefore, the quadratic cost function (21) is concave and the global maximum of the unconstrained problem is obtained for

P = (H + H') (Q^AE^A - Q"E ) (22) 3.3 Constraints

The relative position between the bodies Au(t) is calculated by integrating the difference between (13) and (14); the substitution of Au(t) into (12) provides the expression for the amplitude restriction

Δ"₀ sin («co₀/) ≤ — AU max ,'

which can be written in matrix form as

+ ^Α -X^BY{_∞≤MJ_mai

with

Φ( = [sin (ω₀/), 1 - cos(co₀ , ... , sin (Nco₀t), 1 - cos(Noy )] and W a diagonal matrix W = diag J_ J 1 1_ _1 1_

co₀ ' ω₀ ' 2ω₀ ' 2ω₀ ' ' Nco₀ ' Nco₀

From (1 1 ) and (21 ), it follows that

X^A - X^B = -HP + Q^AE^A - Q^BE^B;

therefore, the amplitude constraint can be expressed as a function of the PTO force and of the excitation forces as: + Φ(ί)(-ΗΡ+ ζ)^ΑΕ^Α -Q^BE^B)W\\ ≤ AU, (23)

4 Optimization algorithm

The algorithm used to solve the constrained optimization problem max J( ) subject to ||Δ«|| < AU_max (24) p

is the penalty method (Nocedal and Wright [2006]); the constrained maximization problem (24) is reformulated as the unconstrained minimization problem

min - J{P) + μ max{0,||A_M| - AU_m , (25) p

where μ > 0 is the penalty parameter. The solution of the problem is sought by solving a sequence of subproblems for increasing values of μ ; the larger μ is the more heavily penalized is the constraint violation in the cost function. If at the step k , where μ = μ^ , the constraint is violated, the subproblem to be solved is:

By the property of the Fourier series (A.4), if there exist a value of P such that max {0,V(2N + l)/2n |Δ«|₂ - AU_max }= 0

that is, if the constraint is satisfied using the 2-norm, then it is also satisfied using the infinity norm:

m^ax |M. -^{A = 0}- The reason for using the 2-norm is that each subproblem min - J(P) + ^(H\₂ - AU_m can be solved analytically and the solution _t ^* is:

P_K' = {H + H^T + 2μ,Η^ΤΨ²ΗΥ (Ι_2Ν + 2μ,Η^ΤΨ²)(0^ΑΕ^Α - Q^BE^B ).

The existence of the solution for the problem (25) is verified by taking the limit

lim P_K' = P: = H- Q^AE^A -Q^BE^B).

The substitution of „^* into (23) results in W AU_max , that is, if the initial condition Au₀ satisfies the constraint then, by increasing μ_¾ , it is possible to find a vector _t ^* for which (23) is satisfied. The limit case is for μ_Α→∞ , where the relative velocity is zero, and the relative position remains constant in the interval / e [Ο, ] .

In practice, the optimization problem is solved by starting with μ, = 0 , which corresponds to the unconstrained problem; if the solution violates the constraint then μ_Α is updated as μ^_+Ι =α μ_¾ with a > 1 , and the new solution P_t ^* ₊₁ is calculated. If the constraint is satisfied the algorithm stops, otherwise the process is repeated until the solution is found. The choice of the parameter a affects both the accuracy and the time required for the computation of the solution. For → 1 the time required for the calculation of the solution increases, because smaller steps are taken, that is the norms (both ||^■ ||₂ and || · \\_∞) decrease slowly. This means that, when increasing μ₄, the first value of _t ^* satisfying the constraint provides a value of the norm ||ΔΗ||_∞ very close to AU_^ . If bigger steps are taken, the first value of P_t ^* may provide ||ΔΜ||_∞ much smaller than AU_rnax, that is restricting the oscillation amplitude more than necessary, thus reducing the produced energy by a significant amount with respect to the optimal value.

The solution provided by the algorithm using the 2-norm is suboptimal when compared to the solution obtained when solving each subproblem using the infinity norm. The main advantage in using the 2-norm is in the time required for the calculation of the solution. 5 Control system implementation

The control system presented in this paper is composed of a feed-forward part and a feedback part; the feed-forward block generates the reference trajectories for the relative velocity, the relative position and the PTO force, that maximize the produced energy while satisfying the amplitude constraint. The feedback controller corrects the PTO force reference signal generated by the feed-forward, in order to minimize the difference between the reference motion and the actual motion of the device. The feedback controller is obtained by solving a continuous time LQ tracking problem as described in Anderson and Moore [1990].

The reference trajectories generated by the feed-forward controller are obtained by solving the optimization problem presented in section 4; the solution of the optimization problem depends on the vectors describing the excitation forces E^A and E^B , which affect both the cost function (21) and the constraint (23). The excitation forces are estimated using the equations of motion (2) and assum ing that measurements of the vertical accelerations and of the relative position are avai lable. Each excitation force is also predicted, using an autoregressive model AR(6), for a prediction horizon of length t .

The reference trajectories are updated periodically with a period At_h ; at every t_h = At_h h , with h = 1 ,2,... , the motion optimization problem and the LQ tracking problem are solved over the time interval T_h = [t_h + t —T,t_h + t ] of length T . The components of the vectors

E^A and E^B are calculated applying the FFT to the time series obtained by the estimation and prediction of the excitation forces over the interval T_h . The length of T_h is such that

T≥ t , that is past values of the excitation forces are also passed to the FFT in order to increase the length of the interval T_h , thus increasing the frequency resolution of the Fourier series. A Tukey window is applied to the signals prior to the FFT to reduce the effect of the spectrum leakage.

6 Simulation results

The control system is simulated using excitation forces calculated from a real sea profile measured by a waverider buoy. The trajectory is updated with a period At_h = 0.5 s, the length of T_h is T = 32 s, the number of frequency components considered is N = 61 and the maximum relative position is set to AU_max = 0.1 m. The simulation results for the relative position Au(t) (Fig. 3) show that when the control is active, from / = 100 s to / = 380 s, the constraint is satisfied. The control algorithm is implemented in MATLAB and the total time required for the computation of the solution of both the constrained motion optimization and the LQ tracking problem is less than 0.3 s on a 2.4GHz dual core personal computer. The time required by the controller for computation can be further reduced by implementing the algorithm on a dedicated hardware.

7 Conclusion

A method for solving of the motion optim ization problem of a self-reacting point- absorber subject to constraint has been presented. The method provides an approximated solution, but allows the constrained optimization problem to be reformulated as a NLP through the discretization of PTO force and of the motion of the bodies. The PTO force and the motion of the device are approximated as linear combinations of basis functions, and the choice of these basis functions determines the properties of the cost function to be optimized and of the constraint. Therefore, the choice of the basis also determines the time required to compute the solution. The approximation by means of the Fourier series has been considered as an example; it provides an orthogonal basis and the resulting cost function is a concave quadratic function . The algorithm used to solve the motion optim ization with ampl itude restriction provides a suboptimal solution, but the convergence is guaranteed and the computational effort is small, making it a candidate for real time implementation.

Different basis functions may be used instead of the Fourier functions, e.g. hat functions or delta functions, and additional constraints can be employed, such as by imposing restrictions on the force applied by the PTO.

Although the device considered is axisymmetric, concentric and oscillating in heave, the method developed can be applied also to a different device oscillating in a different mode.

References Brian D. O. Anderson and John B. Moore, Optimal control: linear quadratic methods. Prentice-Hall, Inc., 1990. ISBNO-l 3-638560-5

. Atkinson and W. Han, Theoretical Numerical Analysis: A Functional Analysis

Framework. Springer, 2005.

J. E. Cuthrell and L. T. Biegler, "On the optimization of differential-algebraic process systems," AIChE J., 33(8): 1257-1270, 1987. ISSN 1547-5905

H. Eidsmoen, "Simulation of a tight-moored amplitude limited heaving-buoy wave- energy converter with phase control," Technical report, NTNU, 1996a.

H. Eidsmoen, "Optimum control of a floating wave-energy converter with restricted

amplitude," Journal of Offshore Mechanics and Arctic Engineering, 1 18(2): 96— 102, 1996b.

D.V. Evans, "Maximum wave-power absorption under motion constraints," Applied Ocean Research, 3(4): 200-203, October 1981 . ISSN 0141- 1 187

A. F. de O. Falcao, P. A. P. Justino, J. C. C. Henriques, and J. M. C. S. Andr. "Reactive versus latching phase control of a two-body heaving wave energy converter," in European Control Conference, 2009.

J. Falnes, "Wave-energy conversion through relative motion between two single-mode oscillating bodies," Journal of Offshore Mechanics and Arctic Engineering, 121(1): 32-38, 1999.

J. Hals, Modelling and phase control of wave-energy converters. PhD thesis, Norwegian University of Science and Technology, NTNU, 2010.

J. Hals, R. Taghipour, and M. Torgeir, "Dynamics of a force-compensated two-body wave energy converter in heave with hydraulic power take-off subject to phase control," in 7th European Wave and Tidal Energy Conference, EWTEC, September 2007.

U. Korde, "Use of oscillation constraints in providing a reaction for deep water floating wave energy devices," International Journal of Offshore and Polar Engineering, 1 1 (2): 1473 - 1484, 2001.

U. A. Korde, "Systems of reactively loaded coupled oscillating bodies in wave energy conversion," Applied Ocean Research, 25(2): 79-91 , Apr. 2003. ISSN 0141 - 1 187 J. Nocedal and S. J. Wright, Numerical Optimization. New York: Springer, 2nd edition,

2006.

T. Ogilvie, "Recent progress toward the understanding and prediction of ship motions," in

Proc. Fifth Symposium on Naval Hydrodynamics, Bergen, Norway, 1964.

A. F. Timan, Theory Of Approximation Of Functions Of A Real Variable. Courier Dover

Publications, 1994.

Appendix A. USEFUL PROPERTIES OF NORMS The 2-norm ||/| of the function / (t) is defined as

||/|₂ = (£l /( I² *)^,e, (A- l) while the infinity-norm ||/||_∞ is defined as

|/|L = sup | /( | . (A.2) ι ,Γ] A general property relating the 2-norm and the infinity-norm of a function / , for which the norm (A.2) exists, is (Atkinson and Han [2005])

\\f\\₂ < fi^₀ ^~ \\fl. (A.3) For a Fourier series with N frequency components f_N , the inequality (Timan [1994])

||_/ <V(2N + l)/2_rt 1_/ (A.4) provides an upper bound for the infinity-norm as a function of the 2-norm.

Appendix B. SOLUTION OF THE EQUATION OF MOTION

For brevity, only the inner products relative to L^A t) will be developed, since the derivations for L⁴ (t) are very similar. Substituting the approximations (5) - (7) into (2) and rearranging, the approximated equation of motion of the body A : L^A ( = Φ (/) - f (t) + S^Au^A (/) (B.1 )

where

r^AA(t) = {m^A +m^A (t) + k^AA(t)*^(i) + B^A +S^A φ,.(τ)Λ

The substitution of the equati the inner product (8) results in:

V j = \,... ,N , where

^ - ]Γ/"( Φ,( Α, gf = ]Γ Φ,( Λ,

0 0

Equation (B.2) can be conveniently rewritten in matrix form as:

G^MX^A +G^ABX^B -T^AP-K^A =0

where g^⁴ are the elements of the matrix G^ , g ^B are the elements of the matrix G^AB ,

— gy are the elements of T^A and -g^A are the elements of K^A .

Equation (10) is then obtained by repeati derivations for the inner products (L^B(t),$_j(t)) and defining the matrix G =

Fig.4 shows a flowchart of operation of a method according to the invention, with

reference to the equations provided above.

The flowchart includes a left-hand offline block 30 and a right-hand online block 32. The offline block 30 represents the programming and calculations performed and the

parameters specific to the particular WEC and PTO for which the control algorithm is

implemented. The online block 32 represents the actual method of the invention as it is operated, i.e. in a loop which continually measures inputs from the WEC system and provides a feedback output to the PTO controller to optimise the operation of the system based on the continually recalculated solution of the optimisation problem.

In block 30, at step 34, an expression is obtained for the velocity and PTO force as a linear combination of basis functions, such as in equations 5-7 or 13-15 (for the Fourier case). This results in step 36 in a set of parameters describing these physical attributes, whose value is not known. In step 38, a matrix version of the equation of motion, e.g. as shown in equation 10 or 18, is generated, and this will, as explained above, depend on a number of physical parameters of the system. In step 40, the cost function and constraint matrices are generated from the equation of motion, e.g. as shown in equations 1 1 and 13 for the general case, or 21 and 23 for the specific Fourier embodiment.

The results of this offline pre-calculation are stored in the computing system used to implement the invention, in the online block 32.

In step 42, the position and acceleration are measured and received as inputs. In step 44, the excitation forces for each body are calculated from these inputs, and in step 46, a prediction of the excitation forces is performed using, for example, an autoregression technique such as AR(6). Step 48 is where the excitation forces are expressed as a combination of basis functions, such as the Fourier series of equations 16 and 17. This results in the parametrisation of the excitation forces in step 50 to provide the parameters for each body. In step 52, these parameters are plugged into the pre-calculated optimisation problem for the cost function and constraint, e.g. as indicated in equation 24 or 25, and as described above in relation to the penalty method, increasing penalty values are applied until a solution is found for the vector P. This provides, in step 54, the sets of velocity parameters for each body and the PTO force parameters, i.e. the missing variables steps 34 and 36. Using a feedback controller in step 56, the PTO force is corrected to match the optimised solution and thereby provide the optimal velocity and position for each body to maximise the energy input. As described previously in the results section of the detailed mathematical treatment, this loop from step 42 to step 56 can be performed in real time to influence and control the operation of the WEC providing maximal energy output while respecting the physical constraints imposed. The loop then continues to repeat in order to maintain the operation of the device continually.

Claims

1 . A computer-implemented method of controlling a power take-off (PTO) of a wave energy converter apparatus having one or more bodies which receive energy from waves in a fluid medium and whose energy is absorbed by said PTO, the method comprising the steps of: receiving an input describing the motion of the or each body making up the system; estimating from said input the excitation force(s) acting on the or each body as a linear combination of basis functions truncated to N terms having coefficients E = [e\, . . . , ^N] ; storing in memory a cost function describing the energy J obtainable from the PTO as a function of a vector P of dimensionality N^p and the coefficients [e_\, ... , ^N] for each body; storing in memory a constraint expression describing a physical constraint on the motion of the or each body, expressed in terms of the vector P and the coefficients [e\, ew] for each body; solving a discrete optimisation problem to maximise the cost function J(P) subject to the constraint expression, thereby providing a solution for the vector P; generating a linear combination of basis functions whose coefficients are the components p_\, ..., p^^p] of the vector P, said linear combination of basis functions providing a force profile for the PTO; and providing as an output to a PTO force controller the calculated PTO force profile.

2. A method as claimed in claim 1 , wherein said cost function is expressed in terms equivalent to a quadratic equation in the vector P.

3. A method as claimed in claim 1 or 2, wherein said cost function is also a function of a set of components [JCI, ...*N] being the coefficients of a linear combination of basis functions representing the velocity of the or each body.

4. A method as claimed in claim 3, wherein said apparatus comprises two bodies A and B, a inndd tthhee vveelloocciittiieess ooff bbooddiieess AA aanndd BB rreessppeectively are approximated as X^A and X^B, where X^A = [x^A, , . . .J ^a _n]^T and X^B = [Λ , ...x^B _N]^T

5. A method as claimed in claim 4, wherein the excitation forces of bodies A and B respectively are approximated as E^A and E^B, where E^A =

...e }.

6. A method as claimed in claim 5, wherein the motion of the system is described by the equation:

X^A

TP + K

X^B where X^A = [x^A · · · ι]^Τ , X^B =

T is an N by N^p matrix, and K is a vector of dimension N; and wherein the cost function J(P) is described by an equation equivalent to: J(P) = P^T (X^A - X^B ) = P^TAT^TAP + AEAP where:

AT = T^I - 1^B;

T⁴ and 7* are obtained by partitioning the matrix G^~'T as [T^T⁶]⁷; ΔΕ = E⁴— E^B; the elements of matrix A are defined as = (Φ, (ΟΦ; (Ο* with { φ i(t), ... , φ Ν(Ϊ)} being the basis of a finite dimensional subspace where the velocities of bodies A and B are approximated as:

v" (0 * v" (0 =∑* , (')

and with { φ ^P|(t), ... , φ ^PNp(t)} being the basis of a finite dimensional subspace where the PTO force is approximated as:

7. A method as claimed in any of claims 3-6, wherein the constraint expression describes a physical constraint on the maximum relative displacement between the two bodies A and B such that the respective vertical positions i (t) and u^B(t) conform to the inequality:

where || . ||_∞ represents the infinity norm.

8. A method as claimed in claim 7, when dependent on claim 6, wherein the stored constraint expression is equivalent to: ||Δ_ΜΟ + Φ(0(ΔΓΡ + Δ£||_¾ < Δ^

where

Δ«₀ = «₀" - u₀ ^B and Φ( = fc⁾, ( .-.9 Af ( ]

and u = 1/(0) ; w₀ ^fi = M^s(0) where

7=1

«^fl( = «„* +∑ ,(/)

7=1

9. A method as claimed in any preceding claim, wherein said basis functions are Fourier series, and wherein said cost function and constraint expression are stored in a form which is specific to the expression of the basis functions as Fourier series.

10. A method as claimed in any preceding claim, wherein the solution of the constrained optimisation problem is performed using a penalty method by solving a sequence of sub-problems with increasing penalty values, until a solution is found which does not violate the constraint.

1 1. A computer program comprising instructions which when executed in a computing system cause the computing system to implement the method of any preceding claim.

12. A computing system programmed with the computer program of claim 1 1 .

13. A wave energy converter comprising one or more bodies which receive energy from waves in a fluid medium, a PTO which applies force to said one or more bodies and thereby absorbs energy from said bodies, one or more sensors for detecting the motion of said one or more bodies, and a computing system as claimed in claim 12, wherein said computing system receives as an input a signal from the one or more sensors and provides as an output a PTO force profile which is implemented by a PTO controller.