Method for establishing the excitation force applied by the swell incident on a movable means of a wave energy system using a model of the drag force

Abstract

The present invention is a method for real-time determination of the forces exerted by incident waves on a mobile part of a wave energy system. Models are constructed of the radiation force exerted on the mobile part and of the drag force exerted on the mobile part and a non-linear model of the wave energy system dynamics. The invention uses only measurements of the float kinematics (position, velocity and possibly acceleration) and of the force applied by a converter machine, which measurements are normally available on a wave energy system since they are used for control and supervision thereof. Determination of the excitation force exerted by incident waves on the mobile part uses the models, the measurements and an unscented Kalman filter.

Claims

1. A computer method for determining an excitation force exerted by incident waves on a mobile part of a wave energy system, the wave energy system converting wave energy to electrical energy through the mobile part cooperating with a conversion machine, which converts the excitation force exerted by the waves on the mobile part into electrical energy, the mobile part moving with respect to the conversion machine under action of the waves, comprising: a) measuring position and velocity of the mobile part by using sensors; b) measuring force exerted by the conversion machine on the mobile part by using sensors; c) constructing a model of radiation force exerted on the mobile part, the model of radiation force relating radiation force to the velocity of the mobile part; d) constructing a model of drag force exerted on the mobile part, the model of drag force relating drag force to velocity of the mobile part; e) constructing a dynamic model of the wave energy system relating the excitation force exerted by the incident waves on the mobile part to position of the mobile part, to velocity of the mobile part, to force exerted by the conversion machine on the mobile part, to the radiation force exerted on the mobile part and to the drag force exerted on the mobile part; and f) determining the excitation force exerted by the incident waves on the mobile part using the dynamic model, the model of radiation force, the model of drag force, the measured position and the measured velocity, and the measured force exerted by the conversion machine on the mobile part, and using an unscented Kalman filter constructed from a random walk model of the excitation force exerted by the incident waves on the mobile part.

2. The computer method as claimed in claim 1, wherein the dynamic model of the wave energy system is constructed by use of an equation: I.sub.eq{umlaut over (δ)}(t)=τ.sub.hd(t)+τ.sub.rad(t)+τ.sub.d(t)+τ.sub.w(t)−u(t), with I.sub.eq being a total moment of inertia of the mobile part, δ(t) being an angle of rotation of the mobile part with respect to an equilibrium position, with {umlaut over (δ)}(t) being an angular acceleration of the mobile part and {dot over (δ)}(t) being an angular velocity of the mobile part, τ.sub.hd(t) being a hydrostatic restoring moment, τ.sub.rad(t) being a radiation moment, τ.sub.d(t) being a drag moment, τ.sub.w(t) being a wave excitation moment and u(t) being a moment exerted by the conversion machine on the mobile part.

3. The computer method as claimed in claim 2, wherein hydrostatic restoring moment τ.sub.hd(t) is determined by use of a formula:
τ.sub.hd(t)=−Kδ(t) where K is a hydrostatic stiffness coefficient.

4. The computer method as claimed in claim 3, wherein the radiation force model is constructed by use of an equation: τ.sub.rad(t)=I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t), with I.sub.∞ being an added moment of inertia at an infinitely high frequency, and τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and {dot over (δ)}(t) being an angular velocity of the mobile part.

5. A method as claimed in claim 4, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

6. A method as claimed in claim 3, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

7. The computer method as claimed in claim 2, wherein the radiation force model is constructed by use of an equation: τ.sub.rad(t)=I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t), with I.sub.∞ being an added moment of inertia at an infinitely high frequency, and τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and {dot over (δ)}(t) being an angular velocity of the mobile part.

8. The computer method as claimed in claim 7, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

9. The computer method as claimed in claim 2, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

10. The computer method as claimed in claim 1, wherein the radiation force model is constructed by use of an equation: τ.sub.rad(t)=−I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t), with I.sub.∞ being an added moment of inertia at an infinitely high frequency, and τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and {dot over (δ)}(t) being an angular velocity of the mobile part.

11. The computer method as claimed in claim 10, wherein the radiation force model is constructed by use of an equation: τ.sub.rad(t)=−I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t), with I.sub.∞ being an added moment of inertia at an infinitely high frequency, and τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and {dot over (δ)}(t) being an angular velocity of the mobile part.

12. The computer method as claimed in claim 10, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

13. A method as claimed in claim 1, wherein the drag force model is constructed by use of an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)| where β is the drag coefficient and {dot over (δ)}(t) is the angular velocity of the mobile part.

14. The computer method as claimed in claim 1, wherein the drag force model F.sub.d is expressed by an equation:
F.sub.d(t)=βż(t)|ż(t)| with β being the drag coefficient and ż being the velocity of the mobile part.

15. The computer method as claimed in claim 1, wherein the wave energy system is controlled depending on the determined excitation force τ.sub.w exerted by the incident waves on the mobile part, by controlling the conversion machine.

16. The computer method as claimed in claim 1, wherein the wave energy system is one of a wave energy system including submerged pressure differential converter device one of a surface floating rotating mass device or a flap device wave energy system.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Other features and advantages of the method according to the invention will be clear from reading the description hereafter, with reference to the accompanying figures wherein:

(2) FIG. 1 illustrates a flop construction wove energy system according to on embodiment of the invention:

(3) FIG. 2 illustrates an estimation of the force exerted by the waves by use of a method according to the prior art for a regular wove;

(4) FIG. 3 illustrates an estimation of the force exerted by the waves by use of a method according to an embodiment of the invention for a regular wave; and

(5) FIG. 4 illustrates an estimation of the force exerted by the waves by use of a method according to an embodiment of the invention for an irregular wove.

DETAILED DESCRIPTION OF THE INVENTION

(6) The present invention relates to a method of determining the excitation force exerted by Incident waves on a mobile port of a wave energy system, which is also referred to as wave excitation force. A wove energy system is a system that converts wave energy to recoverable energy, in particular electrical energy. A wave energy system generally comprises a mobile part, also referred to as flap, pendulum or float, which has an oscillating motion under the action of the waves. The mobile port cooperates with a conversion machine, also referred to as a power take-off (PTO) system, which comprises in most cases an electrical generator coupled to a device allowing adoption of the transmission of the oscillating motion, in order to convert the motion of the mobile part into recoverable energy. In some cases, the conversion machine can act as a motor by generating a force on the mobile port. Indeed, in order to recover power via the conversion machine, a torque or a force is produced opposing the motion of the mobile (generator mode). Furthermore, if the conversion machine allows, power can be supplied to the conversion machine to provide torque or a force for driving the mobile part in order to bring it into resonance with the waves (motor mode).

(7) The method according to the invention is suited for any type of wave energy system with at least one mobile port, as described for example in patent application FR-2,973,448 corresponding to U.S. Pat. No. 9,261,070. The control method according to the invention can also be applied to a wave energy system belonging to the category of wave energy systems with oscillating water columns (OWC). However, the method according to the invention is particularly suited for a non-linear wave energy system, for example of flap construction, of flat submerged pressure differential converter construction or of surface floating rotating mass construction.

(8) FIG. 1 schematically illustrates a non-limitative example of a flap design wave energy system 1. Wave energy system 1 comprises a mobile part 2 in a form of a substantially vertical flap (the wave motion is schematically represented by two curved lines and the direction of the waves is shown by arrow H). Flap 2 has a submerged part and an emerged part. Mobile part 2 cooperates, through a hinge 4, with a converter machine based on an electrical generator 3, which is in this case a rotary generator. Hinge 4 allows a rotary oscillating motion of flop 2 about a horizontal axis with respect to a support 5 that is stationary with respect to the seabed 6.

(9) In the rest of the description, the terms waves, sea waves and wave motion are considered to be equivalent.

(10) Furthermore, in the description, the term force designates a stress or a torque. Similarly, the terms position, velocity and acceleration designate both “linear” and “angular” values. The “linear” values can be associated with the stress and the “angular” values can be associated with the torque. In the rest of the description, only the case relative to torques is illustrated, but the stress-related case can be deduced by transposing the equations to an orthogonal reference frame.

(11) Moreover, for a better understanding, the various models are represented in one dimension. The method according to the invention is suited for models in several dimensions, for systems whose motion has several degrees of freedom.

(12) The method according to the invention comprises the following steps:

(13) 1. Measurement of the mobile part's position and velocity

(14) 2. Measurement of the force exerted by the converter machine (or PTO)

(15) 3. Construction of the radiation force model

(16) 4. Construction of the drag force model

(17) 5. Construction of the dynamic model

(18) 6. Determination of the incident wave excitation force

(19) 7. (optional step) Control of the wave energy system.

(20) Steps 3 and 4 can be carried out in this order, in the reverse order or simultaneously.

(21) The steps of the method can be carried out by a computer (or calculator). The computer can comprise data processing and, advantageously, data storage.

(22) The data processors are configured to implement: measurement of the mobile part's position and velocity. measurement of the force exerted by the converter machine, construction of the radiation force, drag force and dynamic models, and determination of the incident wave excitation force.

(23) 1—Measurement of the Mobile Part's Position and Velocity

(24) This step measures the position and the velocity of the mobile part. The position corresponds to the motion (distance or angle for example) with respect to the position of equilibrium of the mobile part. These measurements can be performed using sensors, generally present on a wave energy system for at least one of control and supervision thereof.

(25) According to an implementation of the invention, this step can also include measuring the acceleration of the mobile part, which can be used to estimate the velocity, or directly in the models used by the method according to the invention. For example, acceleration can be measured using on accelerometer arranged on the mobile part.

(26) 2—Measurement of the Force Exerted by the Converter Machine (PTO)

(27) The force (the stress or possibly the torque) exerted by the PTO converter machine on the mobile part is measured in this step. This measurement can be performed using a sensor, which can be a force sensor or a torque sensor. This type of sensor is often installed or it can be readily installed in wave energy systems, for at least one of control and supervision thereof. Alternatively, the measurement con be replaced by an estimation performed from the force (or torque) setpoint sent to the PTO.

(28) For the wave energy system example illustrated in FIG. 1, a torque sensor can be arranged at hinge 4 (or electrical generator 3).

(29) 3—Construction of the Radiation Force Model

(30) This step constructs a model of the radiation force exerted on the mobile part. According to the linear wave theory (as described for example in the document Folnes J, Kumiawan A. “Fundamental formulae for Wave-Energy Conversion”. R. Soc. open sci. 2: 140305, 2005, http://dx.doi.org/10.1098/rsos.140305), the radiation force results from the oscillation of a submerged body (therefore is depends on the motion of the mobile part), while the excitation force, resulting from the very presence of a body in the water, does not depend on the motion of the submerged body, but is dependent on the incident wave. In the absence of incident wave, the radiation force damps the residual oscillation of the submerged body and eventually stops it. It is important to note that, although the linear theory allows relating the excitation force to the elevation of the incident wave through a linear model (in the frequency or time domain), in practice it cannot be used to calculate the excitation force online, even though it is possible to measure the elevation of the wave at the center of gravity of the mobile part as required by the theory. Indeed, the linear relation between wave elevation and excitation force is non-causal, which means that the excitation force at a given time cannot be calculated without knowing the wave elevation in future times (on the other hand, calculation can be carried out offline, once the wave has passed). In a real-time control context, the excitation force can therefore only be considered as a totally unknown exogenous force acting on the float. On the other hand, still according to the linear wave theory, the radiation force is related to the motion of the float, and more precisely to the acceleration and the velocity thereof, by a causal linear model (in the frequency or time domain). It can therefore be calculated online using the current acceleration measurements and the current and past velocity measurements.

(31) According to an implementation of the invention, the radiation force model τ.sub.rad(t) is constructed by use of an equation:
τ.sub.rad(t)=−I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t)
with:
τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t) the component of the radiation force τ.sub.rad that depends on the (current and past) velocity of the mobile part, which can be referred to as radiation damping;
{umlaut over (δ)} is the angular acceleration of the mobile part;
I.sub.∞ is the added moment of inertia at infinite high frequency, which can be obtained using BEM (Boundary Element Method) calculation codes, such as WAMIT® (WAMIT, USA) or Nemoh® (Ecole Centrale de Nantes, France), from the geometry of the mobile part;
{dot over (δ)} is the angular velocity of the mobile part: and
h is the impulse response relating the velocity of the mobile part to the radiation damping, whose coefficients are obtained from hydrodynamic parameters of the mobile part calculated with the same BEM calculation codes.

(32) The construction of this model allows determination of radiation force at any time, with a limited computation time. Thus, determination of the force exerted by the waves can be determined at any time with a short computation time.

(33) 4—Construction of the Drag Force Model

(34) This step constructs a model of the drag force exerted on the mobile part. The drag force corresponds to the viscous friction on the mobile part of the wave energy system. It is a non-linear model that can depend on the velocity of the mobile port. The drag force (due to viscous friction) is often considered to be negligible for wove energy systems of oscillating point converter type and it is generally excluded from the modelling thereof. This is not the case, however, for flap construction wave energy systems or for machines such as flat submerged pressure differential converters or surface floating rotating mass devices, among others.

(35) According to an embodiment of the invention, for a wave energy system having a rotational motion, the drag force model T.sub.d can be written with an equation:
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)|
with β being the drag coefficient, this coefficient can be determined by experimental tests of the wave energy system or CFD type (Computational Fluid Dynamics) numerical simulations; and
{dot over (δ)} a being the angular velocity of the mobile part.

(36) According to a variant embodiment of the invention, for a wave energy system having a translational motion, the drag force model F.sub.d can be written with an equation:
F.sub.d(t)=βż(t)|ż(t)|
with β being the drag coefficient which can be determined by experimental tests of the wave energy system or CFD type (Computational Fluid Dynamics) numerical simulations, and
ż being the velocity of the mobile part.

(37) The construction of this model enables precise determination of the drag force. Thus, determination of the force exerted by the waves on the mobile part can be precise because it does not disregard the drag force.

(38) 5—Construction of the Dynamic Model of the Wave Energy System

(39) This step constructs a dynamic model of the wave energy system. A dynamic model is understood to be a model that relates the excitation force exerted by the incident waves on the mobile part, the radiation force exerted on the mobile part, the hydrostatic restoring force exerted on the mobile part, the force exerted by the converter machine on the mobile port, to the position and the velocity of the mobile part. This type of model allows obtaining results representative of the behavior of the wave energy system if the motions are not too large.

(40) Advantageously, the dynamic model is obtained by applying the fundamental principle of dynamics (Newton's second law) to the mobile part.

(41) According to an embodiment of the invention wherein the stresses are considered, the dynamic model of the wave energy system can be constructed with an equation:
I.sub.eq{umlaut over (δ)}(t)=τ.sub.hd(t)+τ.sub.rad(t)+τ.sub.d(t)+τ.sub.w(t)−u(t)
with:
I.sub.eq being the moment of inertia of the mobile part,
{umlaut over (δ)} being the angular acceleration of the mobile part,
τ.sub.w being the excitation torque exerted by the incident waves on the mobile part,
τ.sub.rad being the radiation torque exerted on the mobile part,
τ.sub.hd being the hydrostatic restoring torque exerted on the mobile part,
τ.sub.d being the drag torque exerted on the mobile part, and
u being the torque exerted by the converter machine on the mobile part.

(42) This model conveys a rotational motion about a horizontal axis (typical for the wave energy system of FIG. 1). This model is derived from the linear wave theory.

(43) According to a first variant embodiment, the hydrostatic restoring force exerted on the mobile part can be approximated by a linear function of position z defined with respect to the equilibrium position. In this case, the hydrostatic restoring force can be written by a function of the type: τ.sub.hd(t)=−Kδ(t), with δ being the angular position of the mobile part defined with respect to its equilibrium position and K the hydrostatic stiffness coefficient. Thus, the hydrostatic restoring force can be calculated from a simple model if the measurement of position δ is available. This function is particularly well-suited for small displacements δ.

(44) According to an implementation of the invention wherein the stresses are considered, the dynamic model of the wave energy system can be constructed with an equation:
M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)+F.sub.d(t)−F.sub.u(t)
with:
M being the mass of the mobile part,
{umlaut over (z)} being the acceleration of the mobile part,
F.sub.ex being the excitation force exerted by the incident waves on the mobile part,
F.sub.rad being the radiation force exerted on the mobile part,
F.sub.hd being the hydrostatic restoring force exerted on the mobile part,
F.sub.d being the drag force exerted on the mobile part, and
F.sub.u being the force exerted by the converter machine on the mobile part.

(45) This model conveys a vertical translational motion (typical of floats having a heave motion). This model is derived from the linear wave theory.

(46) This model is the “mirror” of the model when the torques are considered; the various terms of the model are similar in nature.

(47) According to an embodiment of the invention, at least one of the hydrostatic restoring force and the hydrostatic restoring torque can be approximated by a linear function or by a piecewise affine function.

(48) 6—Determination of the Excitation Force Exerted by the Incident Waves

(49) In this step, real-time determination of the excitation force exerted by the incident waves on the mobile part is performed by means of: the mobile part's position and velocity (and possibly acceleration) measurements determined in step 1; the measurement of the force exerted by the converter machine PTO on the mobile part determined in step 2; the radiation force model determined in step 3; the drag force model determined in step 4; and the wave energy system dynamic model determined in step 5.

(50) According to the invention, the excitation force exerted by the incident waves on the mobile part is determined using an observer based on an unscented Kalman filter (UKF) constructed from a random walk model of the excitation force exerted by the incident waves on the mobile part. The unscented Kalman filter allows accounting for the non-linearities of the models, in particular the drag force model.

(51) It is noted that a state observer, or state estimator, is, in automation and systems theory, on extension of a model represented as a state-space representation. When the state of the system is not measurable, on observer allowing the state to be reconstructed from a model is constructed.

(52) A UKF filter is based on the “unscented” transformation theory, which allows an estimator to be obtained for a non-linear system without requiring prior linearization for application to the filter. The UKF filter uses a statistical state distribution that is propagated through the non-linear equations. Such a filter has the advantage of providing estimation stability and therefore robustness.

(53) For this embodiment, what is known prior to this step thus is: the mobile part's position δ and velocity δ measurements, the measurement of force u exerted by the conversion machine PTO on the mobile part, the radiation force model τ.sub.rad(t)=−I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t), with τ.sub.r(t)=∫.sub.0.sup.th(t−s){dot over (δ)}(s)ds=h(t)*{dot over (δ)}(t), the drag force model τ.sub.d(t)=δ{dot over (δ)}(t)|{dot over (δ)}(t)|, and the wave energy system dynamic model I.sub.eq{umlaut over (δ)}(t)=τ.sub.hd(t)+τ.sub.rad(t)+τ.sub.d(t)+τ.sub.w(t)−u(t).

(54) In this approach, the problem of estimation of the wave excitation force is transformed into a conventional state estimation problem (that can be solved with an unscented Kalman filter), by expressing the dynamics of the wave excitation force by a random walk model. The main advantage of this method is the consideration of uncertainties allowing accounting for of the measurement noises and the modelled dynamics.

(55) By replacing in the equation that describes the motion of the mobile part
I.sub.eq{umlaut over (δ)}(t)=τ.sub.hd(t)+τ.sub.rad(t)+τ.sub.d(t)+τ.sub.w(t)−u(t)
the expressions for the hydrostatic restoring force (with the linear model), the radiation force model and the drag force model
τ.sub.hd(t)=−Kδ(t)
τ.sub.rad(t)=−I.sub.∞{umlaut over (δ)}(t)−τ.sub.r(t)
τ.sub.d(t)=β{dot over (δ)}(t)|{dot over (δ)}(t)|
the following non-linear model is obtained:
(I.sub.eq+I.sub.∞){umlaut over (δ)}(t)+β{dot over (δ)}(t)|{dot over (δ)}(t)|+Kδ(t)=τ.sub.w(t)+τ.sub.r(t)−u(t).

(56) This equation con be put in state form, by defining

(57) $\{\begin{array}{c}{x}_1\left(t\right)=\delta \left(t\right)\\ x_2\left(t\right)=\delta \left(t\right)\end{array}\hspace{1em}$
which allows the previous model to be written in form of a state-space representation

(58) $\{\begin{array}{c}{\hat{x}}_1\left(t\right)=x_2\left(t\right)\\ {\hat{x}}_1\left(t\right)=\frac{1}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}\left(-{\mathrm{Kx}}_1\left(t\right)-\beta x_2\left(t\right).Math.x_2\left(t\right).Math.-{\tau }_r\left(t\right)+{\tau }_w\left(t\right)-u\left(t\right)\right)\end{array}\hspace{1em}$

(59) This system of equations contains the integral term

(60) ${\tau }_r\left(t\right)={\int }_0^{t}h\left(t-s\right)\stackrel{.}{\delta }\left(s\right)\mathrm{ds}=h\left(t\right)*\stackrel{.}{\delta }\left(t\right)$
that can be considered to be a linear system, with {dot over (δ)}(t) as the input and τ.sub.r(t) as the output. With Prony's method, this system can then be transformed into the equivalent state-space representation:

(61) $\{\begin{array}{c}{\hat{x}}_r\left(t\right)=A_r{x}_r\left(t\right)+B_r\delta \left(t\right)\\ {\tau }_r\left(t\right)=C_r{x}_r\left(t\right)+D_r\delta \left(t\right)\end{array}\mathrm{or}\{\begin{array}{c}{\hat{x}}_r\left(t\right)=A_r{x}_r\left(t\right)+B_r{x}_2\left(t\right)\\ {\tau }_r\left(t\right)=C_r{x}_r\left(t\right)+D_r{x}_2\left(t\right)\end{array}$
where x.sub.r is an internal state (inaccessible) with no particular physical meaning and (A.sub.r, B.sub.r, C.sub.r, D.sub.r) are the matrices of the state-space realization.

(62) By combining the two state-space representations. It is obtained:

(63) $\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_1\left(t\right)=\frac{1}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}\left(-{\mathrm{Kx}}_1\left(t\right)-D_r{x}_2\left(t\right)-C_r{x}_r\left(t\right)-\beta x_2\left(t\right).Math.x_2\left(t\right).Math.+{\tau }_w\left(t\right)-u\left(t\right)\right)\\ {\stackrel{.}{x}}_r\left(t\right)=B_r{x}_2\left(t\right)+A_r{x}_r\left(t\right)\end{array}\hspace{1em}$

(64) or, in an equivalent manner

(65) $\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=A_{c}x\left(t\right)+f_c\left(x\left(t\right)\right)+B_c\left({\tau }_w\left(t\right)-u\left(t\right)\right)\\ y\left(t\right)=C_{c}x\left(t\right)\end{array}\hspace{1em}$

(66) where

(67) $x\left(t\right)={\left[x_1\left(t\right)x_2\left(t\right)x_r^T\left(t\right)\right]}^T$$y=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]$

(68) and

(69) 0$A_c=\left[\begin{array}{ccc}0& 1& 0\\ -\frac{K}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}& -\frac{D_r}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}& -\frac{C_r}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}\\ 0& B_r& A_r\end{array}\right],f_c\left(x\left(t\right)\right)=\left[\begin{array}{c}0\\ -\frac{\beta }{I_{\mathrm{eq}}+I_{\mathrm{oc}}}{x}_2\left(t\right).Math.x_2\left(t\right).Math.\\ 0\end{array}\right]$$B_c=\left[\begin{array}{c}0\\ \frac{1}{I_{\mathrm{eq}}+I_{\mathrm{oc}}}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

(70) Because of the term f.sub.c(x(t)), the system is non-linear. The system has two inputs, τ.sub.w(t) and u(t), and two outputs, x.sub.1(t)=δ(t) and x.sub.2(t)={dot over (δ)}(t). Input τ.sub.w(t) is not measurable and it is unknown. The problem to be solved is to estimate it from the measured quantities u(t) (moment applied by the converter machine PTO), δ(t) and {dot over (δ)}(t) (angular position and velocity of the mobile part).

(71) If the principal motion of the wove energy system were a translational motion and the problem was to estimate the wave excitation force, the some developments would apply by replacing angular position and velocity with positions and velocities, and all the moments with forces.

(72) To carry out the estimation, the above state system is first discretized (because the measurements are sampled and the estimation algorithm is executed by a calculator) using the Euler method, which yields, for a given sampling period T.sub.s:

(73) $\{\begin{array}{c}x\left(k+1\right)=A_{d}x\left(k\right)+f_d\left(x\left(k\right)\right)+B_d{\tau }_w\left(k\right)-B_{d}u\left(k\right)\\ y\left(k\right)=C_{d}x\left(k\right)\end{array}\hspace{1em}$
where A.sub.d=I+T.sub.s, f.sub.d(x(k))=T.sub.sf.sub.c(x(k)), B.sub.d=T.sub.sB.sub.c, C.sub.d=C.sub.c and I is the identity matrix of appropriate dimensions.

(74) To estimate excitation moment τ.sub.w(k), it is considered as a state, by introducing a mathematical model that relates τ.sub.w(k) and τ.sub.w(k+1), in this case:
τ.sub.w(k+1)=τ.sub.w(k)+ϵ.sub.m(k)
where ϵ.sub.m(k) describes the variation of τ.sub.w(k) and is considered to be a random number. In other words, this model implies that, at any time k, the excitation moment strays by a random step (quantity) from its previous value, and that these steps are independently and identically size distributed.

(75) This random walk model of the excitation force is coupled with a more realistic model of the mobile part dynamics:

(76) $\{\begin{array}{c}x\left(k+1\right)=A_{d}x\left(k\right)+f_d\left(x\left(k\right)\right)+B_d{\tau }_w\left(k\right)-B_{d}u\left(k\right)+{ϵ}_x\left(k\right)\\ y\left(k\right)=C_{d}x\left(k\right)+v\left(k\right)\end{array}\hspace{1em}$
where ϵ.sub.x(k) represents the unmodelled dynamics (PTO friction, hydrostatic non-linearity, etc.) and v(k) describes the measurement noise corrupting the float position and velocity measurements.

(77) By combining the random walk model of the excitation force and the non-linear dynamics of the mobile part, the augmented system is obtained:

(78) $\{\begin{array}{c}x\left(k+1\right)=A_{d}x\left(k\right)+f_d\left(x\left(k\right)\right)+B_d{\tau }_w\left(k\right)-B_{d}u\left(k\right)+{ϵ}_x\left(k\right)\\ w\left(k+1\right)={\tau }_w\left(k\right)+{ϵ}_m\left(t\right)\\ y\left(k\right)=C_{d}x\left(t\right)+v\left(k\right)\end{array}\hspace{1em}$
or, in an equivalent manner:

(79) $\{\begin{array}{c}{x}_a\left(k+1\right)=A_a{x}_a\left(k\right)+f_a\left(x_a\left(k\right)\right)+B_{a}u\left(k\right)+ϵ\left(k\right)\\ y\left(k\right)=C_a{x}_a\left(k\right)+v\left(k\right)\end{array}\hspace{1em}$
where

(80) $x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ {\tau }_w\left(k\right)\end{array}\right],f_a\left(x_a\left(k\right)\right)=\left[\begin{array}{c}{f}_d\left(x\left(k\right)\right)\\ 0\end{array}\right],ϵ\left(k\right)=\left[\begin{array}{c}{\mathrm{.Math.}}_x\left(k\right)\\ {\mathrm{.Math.}}_m\left(k\right)\end{array}\right]$
and

(81) $A_a=\left[\begin{array}{cc}{A}_d& B_d\\ 0& 1\end{array}\right],B_a=\left[\begin{array}{c}-B_d\\ 0\end{array}\right],C_a=\left[C_{d}0\right]$

(82) Thus, the problem to be estimated, τ.sub.w(k), becomes a state estimation problem.

(83) One way of estimating the unknown state vector x.sub.a(k), which is by considering information on ϵ(k) and v(k), applies the Kalman filter (KF) algorithm. In this proposed method, the unscented Kalman filter (UKF) is used for handling the non-linearity of the system. The UKF is generally more robust and precise than the extended Kalman filter (EKF), which deals with the non-linearity by linearizing it. Furthermore, in this case, the presence of the drag term, whose derivative is not continuous, makes the EKF inapplicable.

(84) Like the EKF, the UKF performs the estimation in two steps, state prediction and measurement correction, except that these two steps are preceded by a prior step for “sigma points” calculation. The sigma points are a set of samples calculated so as to be able to exactly propagate the mean and variance information in the space of a non-linear function.

(85) According to an implementation of the invention, the following hypotheses can be adopted: the initial state x.sub.a(0) is a random vector of mean m(0)=E[x.sub.a(0)] and of covariance P(0)=E[(x.sub.a(0)−m(0))(x.sub.a(0)−m(0)).sup.T], ϵ(k) and v(k) are Gaussian noises with covariance matrices Q and R respectively,
as well as the following notations: {circumflex over (x)}.sub.a(k|k−1) is the estimation of x.sub.a(k) from measurements up to the time k−1, i.e. y(k−1), y(k−2), . . . and u(k−1), u(k−2), . . . {circumflex over (x)}.sub.a(k|k) is the estimation of x.sub.a(k) from measurements up to the time k, i.e. y(k), y(k−1) . . . and u(k), u(k−1), . . . P.sub.x(k|k−1) is the covariance matrix of x.sub.a(k) from measurements up to the time k−1, i.e. y(k−1), y(k−2), . . . and u(k−1), u(k−2), . . . P.sub.x(k|k) is the covariance matrix of x(k) from measurements up to the time k, i.e. y(k), y(k−1) . . . and u(k), u(k−1), . . .
For this implementation, the three steps in the UKF method can be:

(86) 1. Calculation of the Sigma Points Let:

(87) $W_0^m=\frac{\lambda }{n+\lambda },W_0^c=\frac{\lambda }{n+\lambda }+\left(1-a^2+\gamma \right),W_j^m=W_j^c=\frac{\lambda }{2\left(n+\lambda \right)},j=1,2,\mathrm{.Math.},2n$
where λ=(α.sup.2−1)n is a scaling parameter, n is the dimension of the state x.sub.a(k), α is a parameter that determines the spread of the sigma points around x(k−1|k−1), which is generally assigned a small positive value, 10.sup.−3 for example, γ is a parameter used to incorporate a priori knowledge on the distribution of x: for a Gaussian distribution, γ=2 is optimal.

(88) At the time k−1, the following selection of sigma points (set of points encoding exactly the mean and covariance information is considered):
x.sub.0(k−1)={circumflex over (x)}.sub.a(k−1|k−1),
x.sub.i(k−1)={circumflex over (x)}.sub.a(k−1|k−1)+√{square root over (n+λ)}S.sub.i(k−1), i=1,2, . . . ,n
x.sub.i+n(k−1)={circumflex over (x)}.sub.a(k−1|k−1)−√{square root over (n+λ)}S.sub.i(k−1), i=1,2, . . . ,n
where S.sub.i(k−1) is the i-th column of the matrix square root of P.sub.x(k−1|k−1), that is P.sub.x(k|k−1)=S(k−1).sup.TS(k−1).

(89) 2. Prediction Update

(90) Each sigma point is propagated through the non-linear model representing the evolution of the states:
{circumflex over (x)}.sub.j(k|k−1)=A.sub.ax.sub.j(k−1)+f.sub.a(x.sub.j(k−1))+B.sub.au(k−1), j=0,1, . . . ,2n

(91) The mean and the covariance of {circumflex over (x)}.sub.a(k|k−1), the prediction of x.sub.a(k|k−1) are calculated as

(92) $\hat{x}\left(k|k-1\right)=\underset{j=0}{\overset{2n}{.Math.}}{W}_j^m{\hat{x}}_j\left(k|k-1\right),P_x\left(k|k-1\right)=\underset{j=0}{\overset{2n}{.Math.}}{W}_j^{\sigma }\left({\hat{x}}_j\left(k|k-1\right)-x_a\left(k|k-1\right)\right){\left({\hat{x}}_j\left(k|k-1\right)-x_a\left(k|k-1\right)\right)}^T+Q$

(93) The predicted states {circumflex over (x)}.sub.j(k|k−1) are used in the output state equation, which yields:
ŷ.sub.j(k|k−1)=C.sub.aR.sub.j(k|k−1)

(94) The mean and the covariance of ŷ(k|k−1) are calculated as

(95) $\hat{y}\left(k|k-1\right)=\underset{j=0}{\overset{2n}{.Math.}}{W}_j^m\hat{y}\left(k|k-1\right),P_y\left(k|k-1\right)=\underset{j=0}{\overset{2n}{.Math.}}{W}_j^c\left({\hat{y}}_j\left(k|k-1\right)-y\left(k|k-1\right)\right){\left({\hat{y}}_j\left(k|k-1\right)\right)}^T+R$
while the cross-covariance between {circumflex over (x)}.sub.a(k|k−1) and ŷ(k|k−1) is:

(96) 0$P_{\mathrm{xy}}\left(k.Math.k-1\right)=\underset{j=0}{\overset{2n}{.Math.}}{W}_j^c\left({\hat{x}}_j\left(k.Math.k-1\right)-x_a\left(k.Math.k-1\right)\right){\left({\hat{y}}_j\left(k.Math.k-1\right)-y\left(k.Math.k-1\right)\right)}^T$

(97) 3. Update from the Measurements

(98) As in the linear Kolman filter, the final state estimation is obtained by correcting the prediction with a feedback on the output prediction error (measured):
{circumflex over (x)}(k)={circumflex over (x)}.sub.a(k|k−1)+K(ŷ(k)−ŷ(k|k−1))
where gain K is given by:
K=P.sub.xy(k|k−1)P.sub.y(k|k−1).sup.−1

(99) The a posteriori covariance of the estimation is updated with the formula as follows:
P.sub.x(k|k)=P.sub.x(k|k−1)−KP.sub.y(k|k−1)K.sup.T

(100) The following rules are used to select matrices P.sub.0 (covariance of the initial state) and R: if the initial state {circumflex over (x)}.sub.a(k) at the time k=0 is known, that is m(0)≈x.sub.a(0), then P.sub.0.sup.−1 is large, if there is much noise in the measurements y(k), then R is small.

(101) It is much more complex to select Q. It is generally chosen in diagonal form, as follows:

(102) $Q=\left[\begin{array}{cc}{Q}_x& 0\\ 0& Q_m\end{array}\right]$
with Q.sub.m>>Q.sub.x.

(103) The method based on the UKF and the random walk model of the wave force according to this implementation can be summarized as follows: initializing k=0, state vector {circumflex over (x)}.sub.a(0|0)=m(0) and the state of the covariance matrix P(0|0)=P.sub.0, at any time k: using: the measurements of the mobile part's position and velocity y(k)=[δ(k) {dot over (δ)}(k)] and of the force exerted by the PTO on the mobile part u(k) the results of the estimations of the previous step {circumflex over (x)}.sub.a(k−1|k−1). P(k−1|k−1) parameters Q, R (covariance matrices) determining the excitation force τ.sub.w exerted by the incident waves on the mobile part, denoted by {circumflex over (τ)}.sub.w(k) for this embodiment, by carrying out the following steps: applying the three steps of the unscented Kalman filter algorithm to obtain {circumflex over (x)}.sub.a(k|k), P(k|k) as described above, the complete state being thus estimated with its covariance matrix:
R(k)={circumflex over (x)}.sub.a(k|k−1)+K(ŷ(k)−ŷ(k|k−1))
K=P.sub.xy(k|k−1)P.sub.y(k|k−1).sup.−1
P.sub.x(k|k)=P.sub.x(k|k−1)−KP.sub.y(k|k−1)K.sup.T with x.sub.a(k) being the unknown state vector

(104) $x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ {\tau }_w\left(k\right)\end{array}\right],$
A.sub.a, B.sub.a, C.sub.a, D.sub.a matrices of the state-space realization, P being the covariance matrix of the state vector, Q and R being calibration matrices, and extracting the last component of the estimated state vector, that is the excitation force τ.sub.w exerted by the incident waves on the mobile part, by an equation of the form: ŵ(k)=[0 1] {circumflex over (x)}.sub.a(k|k).

(105) 7—Control of the Wave Energy System

(106) This is an optional step. In this step, the wave energy system is controlled to account for the excitation force exerted by the incident waves. It is thus possible to drive the wave energy system in order to optimize the recovered energy.

(107) This step con control the mobile part of the wave energy system, using for example an electrical, pneumatic or hydraulical conversion machine referred to as PTO (Power Take-Off) system. This PTO system influences the motion of the mobile part and allows the mechanical energy to be transferred to the electrical, pneumatic or hydraulical network. Model predictive control (MPC) is an example of a method of controlling wave energy systems.

COMPARATIVE EXAMPLES

(108) The features and advantages of the method according to the invention will be clear from reading the comparative examples hereafter.

(109) For these examples, a flap wave energy system is used as illustrated in FIG. 1, at small scale for basin testing.

(110) The hydrodynamic parameters of the wave energy system are described in Table 1.

(111) TABLE-US-00001 TABLE 1 Hydrodynamic parameters of the wave energy system Hydrodynamic paramenters Total moment of inertia I.sub.eq + I.sub.∞ 0.2753 kg .Math. m.sup.2 Hydrostatic stiffness coefficient K 118.25 N .Math. m .Math. rad.sup.−1 Drag coefficient β 2 N .Math. m .Math. rad.sup.−1 .Math. s.sup.2

(112) Furthermore, for the state-space representation of the wave energy system, the following matrices are used:

(113) $A_r=\left[\begin{array}{cc}-5\text{,}1899& -36\text{,}7434\\ 1& 0\end{array}\right],B_r=\left[\begin{array}{c}1\\ 0\end{array}\right],C_r=\left[13\text{,}26120\right],D_r=0$

(114) First, the excitation force exerted by the waves on the mobile part is determined for this wave energy system by use of the method according to the prior art described in patent application FR-3,049,989, using a linear Kalman filter after linearization of the drag term. FIG. 2 illustrates the curves of the excitation torque r, in Nm exerted by the waves on the mobile part as a function of time T in s, for a regular wave. Curve REF in solid line corresponds to the reference torque really exerted and curve EST in dotted line corresponds to the torque estimated with this method of the prior art. It is observed that the method according to the prior art with the linear Kalman filter gives imprecise results, with an imprecision that increases with the wave amplitude.

(115) Second, the excitation force exerted by the waves on the mobile part is determined for this wave energy system by use of the method according to the invention (using an unscented Kalman filter). FIG. 3 illustrates the curves of the excitation torque τ.sub.w in Nm exerted by the waves on the mobile part as a function of time T in s, for a regular wave. Curve REF in solid line corresponds to the reference torque really exerted and curve EST in dotted line corresponds to the torque estimated with the method according to the invention. Curves REF and EST are superposed, which shows that the method according to the invention allows precise determination of the excitation torque exerted by the waves on the mobile part.

(116) Thirdly, the second test of FIG. 2 is repeated for irregular waves. FIG. 4 illustrates the curves of the excitation torque r, in Nm exerted by the waves on the mobile port as a function of time T in s, for an irregular wave. Curve REF in solid line corresponds to the reference torque really exerted and curve EST in dotted line corresponds to the torque estimated with the method according to the invention. Curves REF and EST are superposed, which shows that the method according to the invention allows precise determination of the excitation torque exerted by the waves on the mobile part, even for irregular waves.

(117) The method according to the invention is therefore well suited for precisely determining the excitation force exerted by the waves on the mobile part, in particular of a non-linear wave energy system.