METHOD FOR CONTROLLING A WAVE-ENERGY SYSTEM BY DETERMINING THE EXCITATION FORCE APPLIED BY WAVES INCIDENT UPON A MOVING PART OF THE SAID SYSTEM

Abstract

The invention relates to the real-time determination of the forces applied by waves incident upon the moving part (2) of a wave-energy system (1). The method according to the invention is based on the construction of a model of the radiation force applied to the moving part (2) and a model of the dynamics of the wave-energy system (1). The invention uses only measurements of the kinematics of the moving part (2) and the force applied by the conversion machine (3) to the moving part (2).

Claims

1. A method for controlling a wave-energy system by determining the excitation force applied by waves incident upon a moving part of the wave-energy system, the wave-energy system converting the energy of the waves into energy by means of the moving part cooperating with a conversion machine, the moving part performing a movement in relation to the conversion machine as a result of the action of the waves, wherein the following steps are performed: a) measurement of the position and speed of the moving part; b) measurement of the force F.sub.u applied by the conversion machine to the moving part; c) construction of a model of the radiation force F.sub.rad applied to the moving part, the model of the radiation force F.sub.rad relating the radiation force to the speed of the moving part; d) construction of a model of the dynamics of the wave-energy system which relates the excitation force F.sub.ex applied by the wave incident upon the moving part, to the position of the moving part, to the speed of the moving part, to the force F.sub.u applied by the conversion machine to the moving part, and to the radiation force F.sub.rad applied to the moving part; e) determination of the excitation force F.sub.ex applied by the wave incident upon the moving part by means of the dynamics model, the model of the radiation force, the position and speed measured and the measured force F.sub.u applied by the conversion machine to the moving part; and f) controlling of the wave-energy system depending on the determined excitation force F.sub.ex applied by the wave incident upon the moving part.

2. The method as claimed in claim 1, wherein the dynamics model of the wave-energy system is constructed by means of an equation of the type: M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t) where M is the mass of the moving part, {umlaut over (z)} is the acceleration of the moving part, F.sub.hd is the hydrostatic restoring force applied to the moving part.

3. The method as claimed in claim 2, wherein the hydrostatic restoring force F.sub.hd is determined by means of a formula of the type: F.sub.hd(t)=Kz(t) where z is the position of the moving part and K is the hydrostatic stiffness coefficient.

4. The method as claimed in claim 2, wherein the hydrostatic restoring force F.sub.hd is determined by a piecewise affine function of the type: F.sub.hd=F.sub.iz(t)G.sub.i for at least two variation intervals .sub.i of the position z, where F.sub.i and G.sub.i are constant coefficients of the affine function for each interval .sub.i.

5. The method as claimed in claim 1, wherein the model of the radiation force F.sub.rad is constructed by an equation of the type: F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.o.sup.th(t)()d is the radiation damping, {umlaut over (z)} is the acceleration of the moving part, M.sub. is the mass added to the infinitely high frequency, is the speed of the moving part, h is the pulse response which relates the speed of the moving part to the radiation damping.

6. The method as claimed in claim 5, wherein the acceleration of the moving part is measured.

7. The method as claimed in claim 6, wherein the radiation damping F.sub.r is determined by means of an approximated time response obtained from a state realization, of the type: F.sub.r(t)=C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t) where A.sub.r, B.sub.r, C.sub.r, D.sub.r are state realization matrices.

8. The method as claimed in claim 2, wherein the excitation force F.sub.ex applied by the wave incident upon the moving part is determined by means of an equation of the type: {circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.l=0.sup.k-1C.sub.rdA.sub.rd.sup.k-1-lB.sub.rd(l)+D.sub.rd(k)+F.sub.iz(k)+G.sub.i+F.sub.u(k) where k is the time step in the discrete domain, and A.sub.rd, B.sub.rd, C.sub.rd, D.sub.rd are the state realization matrices in the discrete domain.

9. The method as claimed in claim 5, wherein the excitation force F.sub.ex applied by the wave incident upon the moving part is determined by means of a state observer based on a linear Kalman filter constructed from a random-walk model of the excitation force F.sub.ex applied by the said wave incident upon the said moving part.

10. The method as claimed in claim 9, wherein the excitation force F.sub.ex applied by the wave incident upon the moving part is determined by performing the following steps: i) initialization of k=0, the state vector {circumflex over (x)}.sub.a(0|0)=x.sub.a(0) and the state of the covariance matrix, P(0|0)=P.sub.0; ii) at each instant k, acquisition of the position and speed measurements of the moving part y(k)=[z(k) (k)] and the measurement of the force exerted by the conversion machine on the moving part u(k)=F.sub.u(k); and iii) at each instant k, determination of the excitation force applied by the wave incident upon the moving part (2) (k) by means of the following equations:
(k)=[0 1]{circumflex over (x)}.sub.a(k|k) where: $\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a.Math.{\hat{x}}_a\left(k-1|k-1\right)+B_a.Math.F_u\left(k-1\right)\\ P\left(k|k-1\right)=A_a.Math.P\left(k-1|k-1\right).Math.A_a^T+Q\end{array},.Math.\begin{array}{c}K\left(k\right)=P\left(k|k-1\right).Math.{C_a^T\left(C_a.Math.P\left(k|k-1\right).Math.C_a^T+R\right)}^{-1}\\ {\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right).Math.\left(y\left(k\right)-C_a.Math.{\hat{x}}_a\left(k|k-1\right)-D_a.Math.F_u\left(k\right)\right)\\ P\left(k|k\right)=\left(1-K\left(k\right).Math.C_a\left(k\right)\right).Math.P\left(k|k-1\right)\end{array}$ where K(k) is the gain of the observer, x.sub.a(k) is the unknown state vector $x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$ A.sub.a, B.sub.a, C.sub.a, D.sub.a are state realization matrices, P is the covariance matrix of the state vector, and Q and R are calibration matrices.

11. The method as claimed in claim 5, wherein the excitation force FP applied by the wave incident upon the moving part is determined by means of an unknown input state estimator based on a moving-horizon approach.

12. The method as claimed in claim 11, wherein the excitation force F.sub.ex applied by the wave incident upon the moving part is determined by performing the following steps at each instant k: iv) acquisition of the position and speed measurements of the moving part y(k)=[z(k) (k)] and the measurement of the force applied by the conversion machine to the moving part u(k)=F.sub.u(k); and v) determination of the excitation force applied by the wave incident upon the moving part (k|k) by means of the following equations: (k|k)=(N) where (N) is obtained by resolution of a following quadratic programming problem: $\{\begin{array}{c}\underset{}{\min }.Math.\left\{{}^T\left(P+{}_s^T.Math.R.Math..Math.{}_s\right).Math.-2\left[{\hat{w}}_0^T.Math..Math.{\hat{x}}_0^T.Math..Math.V^T\right]\left[\begin{array}{c}S\\ R.Math..Math.{}_s\end{array}\right].Math.\right\}\\ \mathrm{such}.Math..Math.\mathrm{as}.Math..Math.H_{}.Math.g_{}+H_d.Math.U_d\end{array}$ where P, R and S are matrices formed from matrices for weighting of the cost function, is the matrix which represents all the values of the unknown variables in the estimation window of size N, .sub.s is the matrix, calculated from the state-space representation, which corresponds to the effect of these variables on the output values in the same window, H.sub. and H.sub.d are matrices and g.sub. a vector obtained from constraints imposed on the state and the unknown input, {circumflex over (x)}.sub.0={circumflex over (x)}(kN|k1) and .sub.0=(kN|k1) where x is the state vector and w the unknown input of the state estimator, $U_d=\left[\begin{array}{c}{u}_d\left(k-N\right)\\ u_d\left(k-N+1\right)\\ \mathrm{.Math.}\\ u_d\left(k\right)\end{array}\right],$ u.sub.d being the known input of the state estimator with u.sub.d(t)=F.sub.hd(t)F.sub.u(t).

Description

BRIEF DESCRIPTION OF THE FIGURES

[0035] Further characteristic features and advantages of the method according to the invention will become clear from a reading of the description below of non-limiting examples of embodiment, with reference to the attached figures described below.

[0036] FIG. 1 shows wave-energy system according to an embodiment of the invention.

[0037] FIG. 2 shows an approximation of the hydrostatic restoring force as a function of the position according to an example.

[0038] FIG. 3 shows the amplitude of waves for implementation of the illustrative example.

[0039] FIGS. 4a to 4c show the excitation force applied by waves incident upon the moving part as a function of time, for a wave of the type W1 according to FIG. 3 obtained by measurement and respectively by one of the embodiments of the method according to the invention.

[0040] FIGS. 5a to 5c show the excitation force applied by waves incident upon the moving part as a function of time, for a wave of the type W3 according to FIG. 3, obtained by measurement respectively by one of the embodiments of the method according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

[0041] The present invention relates to the determination of the excitation force applied by waves incident upon a moving part of a wave-energy system, called also excitation force of the swell or wave. A wave-energy system is a system which converts the energy of the wave into recoverable energy, in particular electrical energy. A wave-energy system generally comprises a moving part, which is also called pendulum or float, which has an oscillating movement when acted on by the wave. The moving part cooperates with a conversion machine, called also power take-off (PTO) system, which in most cases comprises an electric generator connected to a device allowing the transmission of the oscillating movement to be adapted, so as to convert the movement of the moving part into recoverable energy. In certain cases, the conversion machine may act as a motor generating a force on the moving part. In fact, in order to recover power via a conversion machine, a torque or force which resists the movement of the moving element is produced (generator mode). On the other hand, if the conversion machine allows it, power may be supplied to the conversion machine so as to provide a torque or a force which drives the moving element so as to help it be placed in synchronism with the waves (motor mode). The method according to the invention relates to the control of a wave-energy system depending on the excitation force determined,

[0042] The method according to the invention is suitable for all types of wave-energy system with at least one moving part, for example those described in the patent application FR2973448 (U.S. Pat. No. 9,261,070). The control method according to the invention may also be applied to a wave-energy system belonging to the category of wave-energy systems with oscillating water columns (OWC).

[0043] FIG. 1 shows, in schematic form, a non-limiting example of a wave-energy system 1. The wave-energy system 1 comprises a floating moving part 2 (the wave is represented in schematic form by the hatched area). The moving part 2 cooperates, by means of a lever arm 4, together with a conversion machined based on the electric generator 3 which, in this case, is a linear generator. The lever arm 4 and the electric generator are hingeably mounted relative to an installation 6. The movements of the lever arm are indicated in schematic form by means of arrows.

[0044] In the continuation of the description and the claims, the terms waves, sea tides and swell are considered to be equivalent.

[0045] Moreover, in the description and the claims, the term force indicates a force or a torque. Similarly, the terms position, speed and accelerations all indicate values which are linear or angular, whereby the linear values may be associated with forces, and the angular values may be associated with torques.

[0046] Furthermore, for better comprehension, the different models are shown in a single dimension; however, the method according to the invention is suitable for models with various dimensions, designed for systems which have several degrees of freedom in their movement.

[0047] The method according to the invention comprises the following steps:

[0048] 1) measurement of the position and speed of the moving part;

[0049] 2) measurement of the force applied by the conversion machine (or PTO)

[0050] 3) construction of the radiation force model

[0051] 4) construction of the dynamics model

[0052] 5) determination of the excitation force of the incident wave

[0053] 6) controlling of the wave-energy system

[0054] 1) Measurement of the Position and Speed of the Moving Part:

[0055] During this step, the position and the speed of the moving part are measured. The position corresponds to the movement (for example distance or angle) in relation to the equilibrium position of the moving part. These measurements may be carried out by means of sensors, which are generally provided on a wave-energy system for control and/or management thereof.

[0056] According to a mode of implementation of the invention, during this step, it is also possible to measure or estimate the acceleration of the moving part, since this measurement or estimation may be used in the models implemented by the method according to the invention. For example, the acceleration may be measured by means of an accelerometer mounted on the moving part.

[0057] 2) Measurement of the Force Applied by the Conversion Machine (PTO)

[0058] During this step, the force (or if applicable the torque) applied by the conversion machine (PTO) to the moving part is measured. This measurement may be performed by means of a sensor, which may be a force sensor or a torque sensor. This type of sensor is often installed, or may be easily installed in wave-energy systems, for control and/or management thereof. Alternatively, the measurement may be replaced by an estimation carried out based on a set force (or torque) value sent to the PTO.

[0059] By way of example of the wave-energy system shown in FIG. 1, a torque sensor 5 may be arranged between the lever arm 4 and the electric generator 3.

[0060] 3) Construction of a Radiation Force Model

[0061] During this step, a model of the radiation force applied to the moving part is constructed. According to the linear wave theory (as described for example in the document by Falnes 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 is produced by the oscillation of an immersed body (and therefore depends on the movement of the moving part), while the excitation force, resulting from the presence itself of a body in the water, does not depend on the movement of the immersed body, but on the incident wave. If there is no incident wave, the radiation force dampens the residual oscillation of the immersed body until it causes it to stop. It is important to note that, although the linear theory allows the excitation force to be related to the height of the incident wave by means of a linear model (in the frequency or time domain), in practice it may not be used to calculate the excitation force linearly, even if it was possible to measure the height of the wave at the center of gravity of the float as required by the theory. In fact, the linear relation between height of the wave and excitation force is not random, this meaning that the excitation force cannot be calculated at a given instant without knowing the height of the wave in future instants (the calculation may on the other hand be performed non-linearly once the wave has passed). In the context of real-time control, the excitation force may therefore not be regarded as being 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 movement of the float and more precisely to its acceleration and its speed by a random model which is linear (in the frequency or time domain). It may therefore be calculated linearly using the current acceleration and speed measurements (and past speed measurements).

[0062] According to a mode of implementation of the invention, the said model of the radiation force F.sub.rad is constructed by an equation of the type:

F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t)

[0063] where F.sub.r(t)=.sub.0.sup.th(t)()d is the component of the said radiation force F.sub.rad which depends on the (current and past) speed of the moving part, which may be called radiation damping,

[0064] {umlaut over (z)} is the acceleration of the moving part,

[0065] M.sub. is the mass added to the infinitely high frequency, which may be obtained by means of codes of BEM (Boundary Element Method) calculation codes, such as WAMIT (WAMIT, USA), or Nemoh (Ecole Centrale de Nantes, France), based on the geometry of the moving part,

[0066] is the speed of the moving part,

[0067] h is the pulse response which relates the speed of the moving part to the radiation damping, the coefficients of which are obtained from the hydrodynamic parameters of the moving part calculated by the same BEM calculation codes.

[0068] The construction of this model allows the determination at any instant of the radiation force, with a limited calculation time. Thus the determination of the force applied by the wave may be determined at any instant with a short calculation time.

[0069] 4) Construction of the Dynamics Model of the Wave-Energy System

[0070] During this step, a dynamics model of the wave-energy system is constructed. Dynamics model refers to a model which relates the excitation force applied by waves incident upon the moving part, the radiation force applied to the moving part, the hydrostatic restoring force applied to the moving part and the force applied by the conversion on the moving part, to the position and speed of the moving part. With this type of model it is possible obtain results representing the behavior of the wave-energy system, if the movements are not too big.

[0071] Advantageously, the dynamics model is obtained by applying the fundamental principle of dynamics (Newton's second law) to the moving part.

[0072] According to a mode of implementation of the invention, where the forces are considered, the dynamics model of the wave-energy system may be constructed by an equation of the type:

M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t)

[0073] where M is the mass of the moving part,

[0074] {umlaut over (z)} is the acceleration of the moving part,

[0075] F.sub.ex is the excitation force applied by waves incident upon the moving part,

[0076] F.sub.rad is the radiation force applied to the moving part,

[0077] F.sub.hd is the hydrostatic restoring force applied to the moving part, and

[0078] F.sub.u is the force applied by the conversion machine to the moving part.

[0079] This model translates a vertical translation movement (typically of floats which have a pitching movement). This model is derived from the linear wave theory.

[0080] In accordance with a first variation of embodiment, the hydrostatic restoring force applied to the moving part may be approximated using a linear function of the position z defined in relation to the equilibrium position. In this case, the hydrostatic restoring force may be written by a function of the type: F.sub.hd(t)=Kz(t) where z is the position of the said moving part defined in relation to its equilibrium position and K is the hydrostatic stiffness coefficient. Thus the hydrostatic restoring force may be calculated using a simple model if the measurement of the position z is available. This function is particularly suitable for small displacements z.

[0081] According to a second variation of embodiment, the hydraulic restoring force applied to the moving part may be approximated using a piecewise affine function of the position z. The function may be linear for at least twopreferably threeintervals of the position z. In this case, the hydrostatic restoring force may be written by a function of the type: F.sub.hd=F.sub.iZ(t)G.sub.i for at least two variation intervals .sub.i of the position z, where F.sub.i and G.sub.i are constant coefficients of the affine function for each interval .sub.i. This function is particularly suitable for the larger movements z. FIG. 2 illustrates this approximation. FIG. 2 shows in continuous lines, in a schematic and non-limiting manner, the variations of the hydrostatic restoring force F.sub.hd as a function of the position z. This figure shows in broken lines the linear functions of each piece .sub.1, .sub.2 and .sub.3. It can be noted that this piecewise affine function allows a better approximation of the hydrostatic restoring force over the whole of the variation range of the position z.

[0082] According to a mode of implementation of the invention, where the torques are considered, the dynamics model of the wave-energy system may be constructed by an equation of the type:

J{umlaut over ()}(t)=M.sub.ex(t)+M.sub.hyd(t)+M.sub.rad(t)M.sub.u(t)

[0083] where J is the inertia of the moving part,

[0084] {umlaut over ()} is the angular acceleration of the moving part,

[0085] M.sub.ex is the excitation torque applied by waves incident upon the moving part,

[0086] M.sub.rad is the radiation torque applied to the moving part,

[0087] M.sub.hd is the hydrostatic restoring torque applied to the moving part, and

[0088] M.sub.u is the torque applied by the conversion machine to the moving part.

[0089] This model converts a rotational movement about a horizontal axis (typical of floats which have a pitching movement). This model is derived from the linear wave theory. This model is a mirror of the model where the forces are considered: the various terms of the model are of the same nature.

[0090] In the same way as for the hydrostatic restoring force, the hydrostatic restoring torque may be approximated by a linear function or by a piecewise affine function.

[0091] In the description below, only the example relating to forces will be described, but the example relating to torques may be deduced by means of transposition of the equations in an angular reference system.

[0092] 5) Determination of the Excitation Force Applied by Incident Waves

[0093] During this step, the excitation force applied by waves incident upon the moving part is determined by means of: [0094] measurements of the position and speed (and optionally acceleration) of the moving part, [0095] measurement of the force applied by the conversion machine (PTO) to the moving part, [0096] the model of the radiation force, and [0097] the dynamics model of the wave-energy system.

[0098] According to a first embodiment of the invention (described in detail below), it is possible to determine the excitation force applied by waves incident upon the moving part using the calculation of an approximation of the radiation force. This embodiment has the advantage that a representation of the hydrostatic restoring force by a piecewise affine function may be taken into account.

[0099] According to a second embodiment of the invention (described in detail below), the excitation force applied by waves incident upon the moving part may be determined by means of a state observer based on a linear Kalman filter constructed using a random-walk model of the excitation force applied by waves incident upon the moving part. This second embodiment has the advantage that uncertainty factors are taken into account into account.

[0100] According to a third embodiment of the invention (described in detail below), the excitation force applied by waves incident upon the moving part may be determined by means of an unknown input state estimator based on a moving-horizon approach. This embodiment has the advantage that uncertainty factors are taken into account. Moreover a representation of the hydrostatic restoring force by a piecewise affine function may be taken into account.

[0101] Other embodiments may be envisaged.

First Embodiment

[0102] According to the first embodiment of the invention, it is possible to determine the excitation force applied by waves incident upon the moving part using the calculation of an approximation of the radiation force.

[0103] For this embodiment, the acceleration of the moving part is also measured. Prior to this step, the following are therefore known: [0104] the measurements of the position z, speed and acceleration {umlaut over (z)} of the moving part, [0105] the measurement of the force F.sub.u applied by the conversion machine (PTO) to the moving part, [0106] the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d, and [0107] the dynamics model of the wave-energy system M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t).

[0108] Thus, knowing how to calculate F.sub.r it is therefore possible to obtain F.sub.rad, and finally F.sub.ex(t), the only unknown factor of the dynamics model of the wave-energy system.

[0109] However, it is not easy to perform the calculation of the convolution product within the time domain in the preceding expression for F.sub.r(t), owing to the large number of mathematical operations which must be performed and data to be stored. In order to avoid this, it is possible to regard it as a linear system with F.sub.r(t) at the output and () at the input. The resultant equation in the (Laplace) frequency domain, obtained by applying Prony's method, is:

F.sub.r(s)=W.sub.r(s)(s)

[0110] where W.sub.r(s) is a transfer function. The preceding equation, which was in the Laplace domain, may be expressed in an equivalent state form, for example:

[00005]$\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_r.Math.x_r\left(t\right)+B_r.Math.\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_r.Math.x_r\left(t\right)+D_r.Math.z\left(t\right)\end{array}$

where x.sub.r is an internal state which does not have any particular physical significance and (A.sub.r, B.sub.r, C.sub.r, D.sub.r) are state realization matrices. It is pointed that in system theory (and in automatic control engineering), a state-space representation allows a dynamic system to be modelled in matrix form using state variables. This representation may be linear or non-linear, continuous or discrete. The representation allows determination of the internal state and the outputs of the system at any future instant knowing the state at the initial instant and the behavior of the input variables which influence the system.

[0111] The state form allows the evolution in time of F.sub.r to be expressed as follows:

F.sub.r(t)=C.sub.re.sup.tA.sup.rx.sub.r(0)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)

[0112] but this expression may not be used directly because the initial value of the internal state x.sub.r(0) is unknown. However, it is possible to approximate the time evolution F.sub.r as follows:

F.sub.r(t)=C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)

[0113] since C.sub.re.sup.tA.sup.rx.sub.r(0)0 where t is sufficiently wide, given that the matrix A.sub.r is asymptotically stable for this type of system. This corresponds to ignoring the contribution of the initial state to the time evolution of F.sub.r, knowing that it will recede very rapidly compared to the contribution linked to the speed .

[0114] The first embodiment therefore consists in calculating F.sub.ex from z, , {umlaut over (z)} and F.sub.u (which are known) using an equation with the form:

{circumflex over (F)}.sub.ex(t)=(M+M.sub.){umlaut over (z)}(t)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.r(t)+Kz(t)+F.sub.u(t)

[0115] if a linear model is used for the hydrostatic restoring force, or

{circumflex over (F)}.sub.ex(t)=(M+M.sub.){umlaut over (z)}(t)+C.sub.r.sub.0.sup.te.sup.(t-)A.sup.rB.sub.r()d+D.sub.rz(t)+F.sub.iz(t)+G.sub.i+F.sub.u(t)

[0116] for each interval .sub.i to which z belongs, if a piecewise affine model is used for the hydrostatic restoring force.

[0117] Since in practice the measurements of z, , {umlaut over (z)} and F.sub.u are sampled signals, the algorithm may be developed within the discrete time k, discretizing as follows the state equation introduced beforehand in order to represent the evolution of F.sub.r:

[00006]$\{\begin{array}{c}{x}_r\left(k+1\right)=A_{\mathrm{rd}}.Math.x_r\left(k\right)+B_{\mathrm{rd}}.Math.\stackrel{.}{z}\left(k\right)\\ F_r\left(k\right)=C_{r.Math..Math.d}.Math.x_r\left(k\right)+D_{r.Math..Math.d}.Math.\stackrel{.}{z}\left(k\right)\end{array}$

[0118] F.sub.r may thus be expressed by an equation with the form:

[00007]$F_r\left(k\right)=C_{r.Math..Math.d}.Math.A_{r.Math..Math.d}^k.Math.x_r\left(0\right)+\underset{l=0}{\overset{k-1}{.Math.}}.Math.C_{r.Math..Math.d}.Math.A_{r.Math..Math.d}^{k-1-l}.Math.B_{r.Math..Math.d}.Math.\stackrel{.}{z}\left(l\right)+D_{r.Math..Math.d}.Math.\stackrel{.}{z}\left(k\right)$

[0119] and its approximation expressed by the formula:

[00008]$F_r\left(k\right)=\underset{l=0}{\overset{k-1}{.Math.}}.Math.C_{r.Math..Math.d}.Math.A_{r.Math..Math.d}^{k-1-l}.Math.B_{r.Math..Math.d}.Math.\stackrel{.}{z}\left(l\right)+D_{r.Math..Math.d}.Math.\stackrel{.}{z}\left(k\right)$

[0120] This gives (in the case where the hydrostatic restoring force is considered to be a piecewise affine function) the following expression:

{circumflex over (F)}.sub.ex(k)=(M+M.sub.){umlaut over (z)}(k)+.sub.l=0.sup.k-1C.sub.rdA.sub.rd.sup.k-1-lB.sub.rd(l)+D.sub.rd(k)+F.sub.iz(k)+G.sub.i+F.sub.u(k)

[0121] for each interval .sub.i to which z(k) belongs.

[0122] The original aspect of this embodiment is the approximation of the time response of F.sub.r, the radiation damping (the component of the radiation force dependent on the speed of the float), calculated by means of an equivalent state-space representation and ignoring the contribution of the initial state of this representation.

Second Embodiment

[0123] According to a second embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by means of a state observer based on a linear Kalman filter constructed using a random-walk model of the excitation force applied by waves incident upon the moving part. It is pointed out that a state observer, or state estimator, is in automatic control engineering and system theory an extension of a model represented in the form of a state-space representation. When the system state cannot be measured, an observer allowing the state to be reconstructed from a model is constructed.

[0124] For this embodiment, prior to this step, the following are therefore known: [0125] the measurements of the position z, and speed Z of the moving part, [0126] the measurement of the force F.sub.u applied by the conversion machine (PTO) to the moving part, [0127] the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d, and [0128] the dynamics model of the wave-energy system M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd (t)+F.sub.rad(t)F.sub.u(t).

[0129] In this approach, the problem of estimation of the excitation force of the waves is transformed into a classic state estimation problem (which may be resolved with a linear Kalman filter), where the dynamics of the excitation force of the waves are expressed by a random-walk model. The main advantage of this method is the taking into consideration of uncertainty factors which allows measurement noise and non-modelled dynamics to be taken into account.

[0130] By replacing in the equation which describes the movement of the float

M{umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t)

[0131] the expressions for the hydrostatic restoring force (with the linear model) and the model of the radiation force

F.sub.hd(t)=Kz(t)

F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t)

[0132] the following is obtained

(M+M.sub.){umlaut over (z)}(t)+Kz(t)+F.sub.r(t)=F.sub.ex(t)F.sub.u(t)

[0133] This equation may be expressed in state form, by defining

[00009]$\{\begin{array}{c}{x}_1\left(t\right)=z\left(t\right)\\ x_2\left(t\right)=\stackrel{.}{z}\left(t\right)\end{array}$

[0134] which gives

[00010]$\{\begin{array}{c}.Math.{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ .Math.{\stackrel{.}{x}}_2\left(t\right)=-\frac{K}{M+M_{}}.Math.x_1\left(t\right)-\frac{1}{M+M_{}}.Math.F_r\left(t\right)+\frac{1}{M+M_{}}.Math.\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\end{array}$

[0135] In the same way as for the first embodiment, the dynamics of radiation damping F.sub.r(t) may be described by a state-space representation equivalent to the convolution product of the pulse response of the radiation damping and the speed of the moving part.

F.sub.r(t)=.sub.0.sup.th(t)()d

[0136] namely

[00011]$\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_r.Math.x_r\left(t\right)+B_r.Math.\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_r.Math.x_r\left(t\right)+D_r.Math.\stackrel{.}{z}\left(t\right)\end{array}.Math..Math.\mathrm{or}.Math..Math.\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_r.Math.x_r\left(t\right)+B_r.Math.x_2\left(t\right)\\ F_r\left(t\right)=C_r.Math.x_r\left(t\right)+D_r.Math.x_2\left(t\right)\end{array}$

[0137] where x.sub.r is an internal (non-measurable) state which does not have any particular physical significance, and (A.sub.r, B.sub.r, C.sub.r, D.sub.r) are matrices of the state realization.

[0138] Combining the two state-space representations, the following are obtained:

[00012]$\{\begin{array}{c}.Math.{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ .Math.{\stackrel{.}{x}}_2\left(t\right)=\begin{array}{c}-\frac{K}{M+M_{}}.Math.x_1\left(t\right)-\frac{D_r}{M+M_{}}.Math.x_2\left(t\right)-\frac{C_r}{M+M_{}}.Math.x_r\left(t\right)+\\ \frac{1}{M+M_{}}.Math.\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\end{array}\\ .Math.{\stackrel{.}{x}}_r\left(t\right)=B_r.Math.x_2\left(t\right)+A_r.Math.x_r\left(t\right)\end{array}$

[0139] or, in an equivalent manner

[00013]$\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=A_c.Math.x\left(t\right)+B_c\left(F_{\mathrm{ex}}\left(t\right)-F_u\left(t\right)\right)\\ y\left(t\right)=C_c.Math.x\left(t\right)\end{array}.Math..Math.\mathrm{where}.Math..Math.x\left(t\right)={\left[x_1\left(t\right).Math..Math.x_2\left(t\right).Math..Math.x_r^T\left(t\right)\right]}^T.Math..Math.y=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right].Math..Math.\mathrm{and}.Math..Math.A_c=\left[\begin{array}{ccc}0& 1& 0\\ -\frac{K}{M+M_{}}& -\frac{D_r}{M+M_{}}& -\frac{C_r}{M+M_{}}\\ 0& B_r& A_r\end{array}\right],.Math.B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{}}\\ 0\end{array}\right],C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

[0140] It is necessary to estimate F.sub.ex(t) using y(t) (the vector which contains the position and the speed measurements of the moving part) and F.sub.u(t), the measurement of the force of the PTO.

[0141] In order to do this, firstly the state system above is discretized (it should be remembered that the measurements are sampled), using Tustin's method. This gives

[00014]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d\left(F_{\mathrm{ex}}\left(k\right)-F_u\left(k\right)\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+B_d\left(F_{\mathrm{ex}}\left(k\right)-F_u\left(k\right)\right)\end{array}$

[0142] It should be noted that it is important to perform discretization using Tustin's method because it allows a direct term to appear in the output equation which allows y(k) to be related to F.sub.ex(k), and therefore F.sub.ex(k) to be estimated using y(k). This would not be possible with other discretization methods, which would relate y(k) to F.sub.ex(k1) and not to F.sub.ex(k).

[0143] By defining

w(k)=F.sub.ex(k),u(k)=F.sub.u(k),

[0144] the global model of the float dynamics is obtained:

[00015]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)-B_d.Math.u\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)-D_d.Math.u\left(k\right)\end{array}$

[0145] In order to estimate the excitation force w(k), it is regarded as a state, the evolution of is defined by introducing a mathematical model which relates w(k) to w(k1).

[0146] For this embodiment, a random-walk model is used to describe the evolution of the excitation force:

w(k+1)=w(k)+.sub.m(k)

[0147] where .sub.m(k), which describes the variation of w(k), is regarded as a random number. In other words this model assumes that at each instant k the excitation force deviates by one random step (quantity) from its preceding value and that these steps are distributed independently and identically in terms of size.

[0148] This random-walk model of the excitation force is then combined with a dynamics model of the float:

[00016]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)-B_d.Math.u\left(k\right)+{}_x\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)-D_d.Math.u\left(k\right)+\left(k\right)\end{array}$

[0149] which is more resistant to modelling errors, owing to the introduction of .sub.x(k) which represents the unmodelled dynamics (friction of the PTO, hydrostatic non-linearity, . . . ) and (k) which describes the measurement noise affecting the position and speed of the float.

[0150] By combining the random-walk model of the excitation force of the incident waves and the dynamics of the float. the increased state system is obtained:

[00017]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)-B_d.Math.u\left(k\right)+{}_x\left(k\right)\\ w\left(k+1\right)=w\left(k\right)+{}_m\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)-D_d.Math.u\left(k\right)+\left(k\right)\end{array}$

[0151] or, in an equivalent manner

[00018]$\{\begin{array}{c}{x}_a\left(k+1\right)=A_a.Math.x_a\left(k\right)+B_a.Math.u\left(k\right)+\left(k\right)\\ y\left(k\right)=C_a.Math.x_a\left(k\right)+D_a.Math.u\left(k\right)+\left(k\right)\end{array}.Math..Math.\mathrm{where}.Math..Math.x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],\left(k\right)=\left[\begin{array}{c}{}_x\left(k\right)\\ {}_m\left(k\right)\end{array}\right].Math..Math.\mathrm{and}.Math..Math.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[\begin{array}{cc}{C}_d& D_d\end{array}\right],D_a=-D_d$

[0152] Thus the problem of estimating w(k) becomes a state estimation problem.

[0153] One way of estimating the unknown state vector x.sub.a(k), taking into account information about (k) and (k), is to apply the linear Kalman filter algorithm, as shown below. An alternative may be to use the moving-horizon estimation approach (the same approach used in the third embodiment to estimate the excitation force, viewed this time as an unknown input).

[0154] The linear Kalman filter (KF) algorithm provides a solution for the following optimization problem:

[00019]$\underset{x_a\left(k\right)}{\min }.Math.=\left\{{\left(x_a\left(0\right)-{\stackrel{\_}{x}}_a\left(0\right)\right)}^T.Math.P_0^{-1}\left(x_a\left(0\right)-{\stackrel{\_}{x}}_a\left(0\right)\right)+\underset{j=1}{\overset{k}{.Math.}}.Math.{\left(j-1\right)}^T.Math.Q^{-1}.Math.\left(j-1\right)+{\left(j\right)}^T.Math.R^{-1}.Math.\left(j\right)\right\}$

[0155] where P.sub.0, Q and R are calibration matrices of suitable dimensions, {circumflex over (x)}.sub.a(0) is the mean value of the initial state x.sub.a(0).

[0156] Considering x.sub.a(0) to be a random vector, not related to (k) and (k), of mean {circumflex over (x)}.sub.a(0) with covariance matrix P.sub.0, and (k) and (k) as two random processes of the mean-zero white noise type with covariance matrices Q and R respectively, the problem of optimization may be solved using the linear Kalman filter algorithm.

[0157] At each instant k, the linear Kalman filter algorithm calculates the solution to this problem by means of two steps.

[0158] The first step is the temporal updating of the estimations:

[00020]$\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a.Math.{\hat{x}}_a\left(k-1|k-1\right)+B_a.Math.F_u\left(k-1\right)\\ P\left(k|k-1\right)=A_a.Math.P\left(k-1|k-1\right).Math.A_a^T+Q\end{array}$

[0159] where {circumflex over (x)}.sub.a(k|k1) and P(k|k1) are respectively the estimation of x.sub.a(k) and its covariance matrix obtained using the measurements from the instant k1 (namely y(k1), y(k2), . . . and F.sub.u(k1), F.sub.u(k2), . . . ), and {circumflex over (x)}.sub.a(k|k) and P(k|k) are the estimation of x.sub.a(k) and its covariance matrix obtained using the measurements from the instant k (namely y(k), y(k1), . . . and F.sub.u(k), F.sub.u(k1), . . . ).

[0160] The second step is the updating of the measurements:

[00021]$\{\begin{array}{c}K\left(k\right)=P\left(k|k-1\right).Math.{C_a^T\left(C_a.Math.P\left(k|k-1\right).Math.C_a^T+R\right)}^{-1}\\ {\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right).Math.\left(y\left(k\right)-C_a.Math.{\hat{x}}_a\left(k|k-1\right)-D_a.Math.M_u\left(k\right)\right)\\ P\left(k|k\right)=\left(I-K\left(k\right).Math.C_a\left(k\right)\right).Math.P\left(k|k-1\right)\end{array}$

[0161] with K(k) the gain of the observer, which is a parameter vector which weights the errors between the estimations and the measurements of the outputs of the system (position and speed). In the case of the Kalman filter, it is not a calibration parameter, but it is obtained from the estimation of the covariance matrix.

[0162] To summarize, for the second embodiment, the following steps may be performed: [0163] Initialization of k=0, {circumflex over (x)}.sub.a(0|0)={circumflex over (x)}.sub.a(0), P(0|0)=P.sub.0 [0164] At each instant k: [0165] the following are used: [0166] the measurements of the position and speed of the moving part y(k)=[z(k) (k)] and the force applied by the PTO to the moving part u(k)=F.sub.u (k) [0167] the results of the estimations of the preceding step {circumflex over (x)}.sub.a(k1|k1), P(k1|k1) [0168] the parameters Q, R (covariance matrices) [0169] determination of the excitation force F.sub.ex applied by the waves incident upon the moving part, known for this embodiment (k), by performing the following steps: [0170] the two steps of the linear Kalman filter algorithm are applied to obtain {circumflex over (x)}.sub.a(k|k), P(k|k), thus the complete state with its covariance matrix are estimated:

[00022]$\{\begin{array}{c}{\hat{x}}_a\left(k|k-1\right)=A_a.Math.{\hat{x}}_a\left(k-1|k-1\right)+B_a.Math.F_u\left(k-1\right)\\ P\left(k|k-1\right)=A_a.Math.P\left(k-1|k-1\right).Math.A_a^T+Q\end{array},.Math.\begin{array}{c}K\left(k\right)=P\left(k|k-1\right).Math.{C_a^T\left(C_a.Math.P\left(k|k-1\right).Math.C_a^T+R\right)}^{-1}\\ {\hat{x}}_a\left(k|k\right)={\hat{x}}_a\left(k|k-1\right)+K\left(k\right).Math.\left(y\left(k\right)-C_a.Math.{\hat{x}}_a\left(k|k-1\right)-D_a.Math.F_u\left(k\right)\right)\\ P\left(k|k\right)=\left(I-K\left(k\right).Math.C_a\left(k\right)\right).Math.P\left(k|k-1\right)\end{array}$ [0171] where x.sub.a(k) is the unknown state vector

[00023]$x_a\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ w\left(k\right)\end{array}\right],$ A.sub.a, B.sub.a, C.sub.a, D.sub.a are state realization matrices, P is the covariance matrix of the state vector, and Q and R are calibration matrices. [0172] and extraction of the last component of the last component of the estimated state vector, namely the excitation force F.sub.ex applied by the waves incident upon the moving part, by means of an equation with the form: (k)=[0 1] {circumflex over (x)}.sub.a(k|k)

Third Embodiment

[0173] According to the third embodiment of the invention, the excitation force applied by waves incident upon the moving part may be determined by means of an unknown input state estimator based on a moving-horizon approach. The moving-horizon estimation is an algorithm which determines recursively the state and the unknown input using data belong to a finite time window.

[0174] For this embodiment, prior to this step, the following are therefore known: [0175] the measurements of the position z, and speed of the moving part, [0176] the measurement of the force F.sub.u applied by the conversion machine (PTO) to the moving part, [0177] the model of the radiation force F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t) where F.sub.r(t)=.sub.0.sup.th(t)()d, and [0178] the dynamics model of the wave-energy system {umlaut over (z)}(t)=F.sub.ex(t)+F.sub.hd(t)+F.sub.rad(t)F.sub.u(t).

[0179] The second embodiment is completely linear and does not allow a non-linear (more precise) description of the hydrostatic restoring force to be taken into account directly.

[0180] This third embodiment allows this limitation to be overcome. Unlike the second embodiment, the excitation force is no longer regarded as a global state model which represents the dynamics of the float, the evolution of which is described by a random-walk model, but is regarded as an unknown input of a (slightly different) global model representing the dynamics of the float.

[0181] By combining the equation which describes the movement of the float according to the linear wave theory

M{umlaut over (z)}(t)=F.sub.ex(t)+F(t)+F.sub.rad(t)F.sub.u(t)

[0182] and the expression for the radiation force

F.sub.rad(t)=M.sub.{umlaut over (z)}(t)F.sub.r(t)

[0183] the following is obtained

(M+M.sub.){umlaut over (z)}(t)+F.sub.r(t)=w(t)+u.sub.d(t)

[0184] where w(t)=F.sub.ex(t) and u.sub.d(t)=F.sub.hd(t)F.sub.u(t). u.sub.d(t) is regarded as a control input of the system described, known because F.sub.u(t) is measurable and F.sub.hd(t) may be calculated from a from a (linear or non-linear) model of the hydrostatic restoring force, while w(t) is an unknown input (disturbance).

[0185] This equation may be expressed in state form, by defining

[00024]$\{\begin{array}{c}{x}_1\left(t\right)=z\left(t\right)\\ x_2\left(t\right)=\stackrel{.}{z}\left(t\right)\end{array}$

[0186] which gives

[00025]$\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{1}{M+M_{}}.Math.F_r\left(t\right)+\frac{1}{M+M_{}}.Math.\left(w\left(t\right)-u_d\left(t\right)\right)\end{array}$

[0187] In the same way as for the second embodiment, using

F.sub.r(t)=.sub.0.sup.th(t)()d

[0188] the radiation damping F.sub.r(t) may be expressed as the output of a linear system in the form of a state-space representation, the input of which is x.sub.2(t)=(t)

[0189] either

[00026]$\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_r.Math.x_r\left(t\right)+B_r.Math.\stackrel{.}{z}\left(t\right)\\ F_r\left(t\right)=C_r.Math.x_r\left(t\right)+D_r.Math.\stackrel{.}{z}\left(t\right)\end{array}.Math..Math.\mathrm{or}.Math..Math.\{\begin{array}{c}{\stackrel{.}{x}}_r\left(t\right)=A_r.Math.x_r\left(t\right)+B_r.Math.x_2\left(t\right)\\ F_r\left(t\right)=C_r.Math.x_r\left(t\right)+D_r.Math.x_2\left(t\right)\end{array}$

[0190] where x.sub.r is a (non-measurable) internal state not having any particular physical significance and (A.sub.r, B.sub.r, C.sub.r, D.sub.r) are the state realization matrices By combining the two state-space representations, the following is obtained:

[00027]$\{\begin{array}{c}{\stackrel{.}{x}}_1\left(t\right)=x_2\left(t\right)\\ {\stackrel{.}{x}}_2\left(t\right)=-\frac{D_r}{M+M_{}}.Math.x_2\left(t\right)-\frac{C_r}{M+M_{}}.Math.x_r\left(t\right)+\frac{1}{M+M_{}}.Math.\left(w\left(t\right)-u_d\left(t\right)\right)\\ {\stackrel{.}{x}}_r\left(t\right)=B_r.Math.x_2\left(t\right)+A_r.Math.x_r\left(t\right)\end{array}$

[0191] or, in an equivalent manner

[00028]$\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=A_c.Math.x\left(t\right)+B_c\left(w\left(t\right)+u_d\left(t\right)\right)\\ y\left(t\right)=C_c.Math.x\left(t\right)\end{array}.Math..Math.\mathrm{where}.Math..Math.x\left(t\right)={\left[\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_r^T\left(t\right)\end{array}\right]}^T.Math..Math.y=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right].Math..Math.\mathrm{and}.Math..Math.A_c=\left[\begin{array}{ccc}0& 1& 0\\ 0& -\frac{D_r}{M+M_{}}& -\frac{C_r}{M+M_{}}\\ 0& B_r& A_r\end{array}\right],.Math.B_c=\left[\begin{array}{c}0\\ \frac{1}{M+M_{}}\\ 0\end{array}\right],.Math.C_c=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

[0192] It is attempted to estimate w(t) using y(t) (the vector which contains the measurements of the position and the speed of the float) and u.sub.d(t). It is pointed out that u.sub.d(t)=F.sub.hd(t)F.sub.u(t) and that it is possible to use, for example, the piecewise affine approximation in order to calculate the hydrostatic restoring force:

F.sub.hd=F.sub.iz(t)G.sub.i for each interval .sub.i, to which z belongs.

[0193] To do this, it is first possible to discretize the abovementioned state system (it should be remembered the measurements are sampled), using Tustin's method. This gives:

[00029]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)+B_d.Math.u_d\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)+D_d.Math.u_d\left(k\right)\end{array}$

[0194] To make the model more realistic by adding uncertainty factors (noise) to the state and the output:

[00030]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)+B_d.Math.u_d\left(k\right)+\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)+D_d.Math.u_d\left(k\right)+\left(k\right)\end{array}$

[0195] where (k) and (k) are mean-zero white noise.

[0196] In order to estimate the unknown input w(k) (and secondarily the state x(k)) the moving-horizon estimation (MHE) method is used. The problem consists in estimating, at each time instant kN, the state vector x(k) and the unknown input w(k)=F.sub.ex(k), on the basis of: [0197] a priori estimations: {circumflex over (x)}(kN|k1), (kN|k1) (the estimations of x(kN) and w(kN) respectively) using the measurements from the time k1; [0198] information about the outputs: y(kN), y(kN+1), . . . , y(k); and [0199] information about the known input: u.sub.d(kN), u.sub.d(kN+1), . . . , u.sub.d(k)

[0200] Namely:

{circumflex over (x)}.sub.0={circumflex over (x)}(kN|k1)

.sub.0=(kN|k1)

and:

x.sub.j={circumflex over (x)}(kN+j|k)

w.sub.j=(kN+j|k)

.sub.j={circumflex over ()}(kN+j|k)

[0201] for j=0, 1, . . . , N.

[0202] The moving-horizon estimation problem is formulated as:

[00031]$\underset{x_0,w_0,w_1,\mathrm{.Math.}.Math.,{}_1,\mathrm{.Math.}.Math.,{}_N}{\min }.Math.J\left(x_0,w_0,w_1,\mathrm{.Math.}.Math.,{}_1,\mathrm{.Math.}.Math.,{}_N\right)$

[0203] with the following objective function (cost):

[00032]$J=(.Math.)={.Math.x_0-{\hat{x}}_0.Math.}_{P_0^{-1}}^2+{.Math.w_0-{\hat{w}}_0.Math.}_{Q_0^{-1}}^2+\underset{j=0}{\overset{N-1}{.Math.}}.Math.{.Math.{}_j.Math.}_{Q^{-1}}^2++.Math.\underset{j=0}{\overset{N}{.Math.}}.Math.{.Math.y\left(k-N+j\right)-\hat{y}\left(k-N+j\right).Math.}_{R^{-1}}^2+\underset{j=0}{\overset{N-1}{.Math.}}.Math.{.Math.w_j-w_{j-1}.Math.}^2$

[0204] where

(kN+j)=C.sub.dx.sub.j+D.sub.dw.sub.j+D.sub.du.sub.d(kN+j)

[0205] The notation v.sub.M.sup.2 indicates the product v.sup.TMv for a vector of reals v and a symmetrical matrix of reals M.

[0206] There are five terms in the cost function: [0207] The first term

[00033]${.Math.x_0-{\hat{x}}_0.Math.}_{P_0^{-1}}^2$

and a second term

[00034]${.Math.w_0-{\hat{w}}_0.Math.}_{Q_0^{-1}}^2$

added together form the arrival cost which summarizes information prior to the time kN and forms part of the data of the estimation problem. [0208] The third term and the fourth term penalize the method and measurement noise, respectively, with the weighting matrices Q and R. [0209] The last term is a regularization term which has the function of taking into account the fact that the excitation force of the waves is a regular (smooth) signal.

[0210] This latter term is important for determining the force applied by waves in an optimum manner. Another way of taking into account the smooth nature of the excitation force of waves would consist in imposing variation constraints thereon.

[0211] Using the classic approach of a moving-horizon estimation and control, the minimization of the cost function may be transformed into a quadratic programming QP problem (a quadratic programming problem is an optimization problem in which a quadratic function is minimized/maximized on a convex polyhedron) by propagating:

[00035]$\{\begin{array}{c}x\left(k+1\right)=A_d.Math.x\left(k\right)+B_d.Math.w\left(k\right)+B_d.Math.u_d\left(k\right)+\left(k\right)\\ y\left(k\right)=C_d.Math.x\left(k\right)+D_d.Math.w\left(k\right)+D_d.Math.u_d\left(k\right)+\left(k\right)\end{array}$ [0212] over time. This gives:

[00036]$.Math.x_1=A_d.Math.x_0+B_d.Math.w_0+B_d.Math.u_d\left(k-N\right)+{}_0$$\begin{array}{c}.Math.x_2=.Math.A_d.Math.x_1+B_d.Math.w_1+B_d.Math.u_d\left(k-N+1\right)+{}_1\\ =.Math.A_d^2.Math.x_0+A_d.Math.B_d.Math.w_0+B_d.Math.w_1+A_d.Math.B_d.Math.u_d\left(k-N\right)+\\ .Math.B_d.Math.u_d\left(k-N+1\right)+A_d.Math.{}_0+{}_1\end{array}$$.Math.\mathrm{.Math.}$$x_n=A_d^N.Math.x_0+A_d^{N-1}.Math.B_d.Math.w_0+\mathrm{.Math.}+B_d.Math.w_{N-1}+A_d^{N-1}.Math.B_d.Math.u_d\left(k-N\right)+\mathrm{.Math.}+B_d.Math.u_d\left(k-1\right)+A_d^{N-1}.Math.{}_0+\mathrm{.Math.}+{}_1$ [0213] or, in compact form,

[00037]$X={}_x.Math.x_0+{}_x.Math.W+{}_x.Math.U_d+{}_x.Math.E$$\mathrm{where}$$X=\left[\begin{array}{c}{x}_1\\ x_2\\ \mathrm{.Math.}\\ A_N\end{array}\right],W=\left[\begin{array}{c}{w}_1\\ w_2\\ \mathrm{.Math.}\\ W_N\end{array}\right],U_d=\left[\begin{array}{c}{u}_d\left(k-N\right)\\ u_d\left(k-N+1\right)\\ \mathrm{.Math.}\\ u_d\left(k\right)\end{array}\right],E=\left[\begin{array}{c}{}_1\\ {}_2\\ \mathrm{.Math.}\\ {}_N\end{array}\right]$$\mathrm{and}$${}_x=\left[\begin{array}{c}{A}_d\\ A_d^2\\ \mathrm{.Math.}\\ A_d^N\end{array}\right],{}_x=\left[\begin{array}{ccccc}{B}_d& 0& \mathrm{.Math.}& 0& 0\\ A_d.Math.B_d& B_d& \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ A_d^{N-1}.Math.B_d& A_d^{N-2}.Math.B_d& \mathrm{.Math.}& B_d& 0\end{array}\right],.Math.{}_x=\left[\begin{array}{cccc}I& 0& \mathrm{.Math.}& 0\\ A_d& I& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ A_d^{N-1}& A_d^{N-2}& \mathrm{.Math.}& I\end{array}\right]$

[0214] Applying the same approach to the output equation and introducing the notations

[00038]$Y={\left[{y\left(k-N+1\right)}^T.Math..Math.{y\left(k-N+2\right)}^T.Math..Math.{y\left(k\right)}^T\right]}^T,.Math.{}_y=\left[\begin{array}{c}{C}_d.Math.A_d\\ C_d.Math.A_d^2\\ \mathrm{.Math.}\\ C_d.Math.A_d^N\end{array}\right],.Math.{}_y=\left[\begin{array}{ccccc}{C}_d.Math.B_d& D_d& \mathrm{.Math.}& 0& 0\\ C_d.Math.A_d.Math.B_d& C_d.Math.B_d& \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ C_d.Math.A_d^{N-1}.Math.B_d& C_d.Math.A_d^{N-2}.Math.B_d& \mathrm{.Math.}& C_d.Math.B_d& D_d\end{array}\right],.Math.{}_y=\left[\begin{array}{cccc}{C}_d& 0& \mathrm{.Math.}& 0\\ C_d.Math.A_d& C_d& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ C_d.Math.A_d^{N-1}& C_d.Math.A_d^{N-2}& \mathrm{.Math.}& C_d\end{array}\right]$

[0215] results in a QP problem with the form:

[00039]$\{\begin{array}{c}\underset{}{\min }.Math.\left\{{}^T\left(P+{}_s^T.Math.R.Math..Math.{}_s\right).Math.-2\left[{\hat{w}}_0^T.Math..Math.{\hat{x}}_0^T.Math..Math.V^T\right]\left[\begin{array}{c}S\\ R.Math..Math.{}_s\end{array}\right].Math.\right\}\\ \mathrm{such}.Math..Math.\mathrm{as}.Math..Math.H_{}.Math.g_{}+H_d.Math.U_d\end{array}.Math..Math.\mathrm{where}.Math..Math.\{\begin{array}{c}=\left[W^T.Math..Math.x_0^T.Math..Math.E^T\right]\\ V=Y-{}_y.Math.U_d\\ {}_s=\left[{}_y.Math..Math.{}_y.Math..Math.{}_y\right]\end{array}$

[0216] where is a matrix which represents all the unknown variables to be estimated (the values of the unknown input W.sup.T and the uncertainty about the state E.sup.T for all the N steps of the estimation window considered, and the initial state x.sub.0.sup.T at the time kN), .sub.s is a matrix which includes the effects of these same unknown variables on the values output in the estimation window Y, .sub.y, is a matrix which corresponds to the effects of the known input U.sub.d on Y, .sub.y .sub.y and .sub.s calculated based on the state-space representation defined above, and P, R and S are matrices formed from weighting matrices of the cost function Q, R, Q.sub.0.sup.1, R.sub.0.sup.1 and the regularization parameter as follows:

[00040]$Q=\left[\begin{array}{cccc}{Q}^{-1}& 0& \mathrm{.Math.}& 0\\ 0& Q^{-1}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& Q^{-1}\end{array}\right],R=\left[\begin{array}{cccc}{R}^{-1}& 0& \mathrm{.Math.}& 0\\ 0& R^{-1}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& R^{-1}\end{array}\right],.Math.{}_1=\left[\begin{array}{ccccccc}{Q}_0^{-1}-& & 0& 0& \mathrm{.Math.}& 0& 0\\ 0& -& & 0& \mathrm{.Math.}& 0& 0\\ 0& 0& -& & \mathrm{.Math.}& 0& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}& \mathrm{.Math.}\\ 0& 0& 0& 0& \mathrm{.Math.}& & 0\\ 0& 0& 0& 0& \mathrm{.Math.}& -& \end{array}\right],.Math.P=\left[\begin{array}{ccc}{}_1& 0& 0\\ 0& P_0^{-1}& 0\\ 0& 0& Q\end{array}\right],S=\left[\begin{array}{cccc}{Q}_0^{-1}& 0& \mathrm{.Math.}& 0\\ 0& P_0^{-1}& \mathrm{.Math.}& 0\end{array}\right]$

[0217] To conclude

H.sub.g.sub.+H.sub.dU.sub.d

[0218] is a linear inequality which, with the aid of the matrices H.sub. and H.sub.d and the vector g.sub., retranscribes in compact form polytopic constraints on the unknown input and state:

[00041]$\{\begin{array}{c}w{P}_w,P_w=\left\{w.Math.\text{:}.Math..Math.H_w{g}_w\right\},g_w>0\\ x{P}_x,P_x=\left\{x.Math.\text{:}.Math..Math.H_x{g}_x\right\},g_x>0\end{array}$

[0219] which in turn are defined by the matrices H.sub.w and H.sub.x and the vectors g.sub.w and g.sub.x. With this inequality it is possible to take into account information about the minimum and maximum amplitudes of the excitation force of the incident waves, and the position and speed of the float.

[0220] The moving-horizon estimation problem is a constrained quadratic programming problem which may be resolved with solvers such as FORCES (Embotech, Suisse) or CVXGEN (CVXGEN LLC, USA) for example.

[0221] Let it be assumed that * is the solution of the dynamic programming problem obtained with one of these solvers. The estimation of the excitation forces, using measurements at the time k, is given by:

(kN+j|k)=*(j),i=1,2, . . . ,N

[0222] To summarize, for the third embodiment, the steps described below may be performed.

[0223] The components of the vectors of the past estimations and measurements corresponding to the N1 time step preceding the present instant k=0 are initialized to zero.

[0224] At each instant k, [0225] 1. The following are acquired: [0226] the measurements of the position and the speed of the moving part y(k)=[z(k) (k)] and of the force applied by the PTO to the moving part u(k)=F.sub.u (k), [0227] the results of the estimations of the preceding time step {circumflex over (x)}.sub.d(kN|k1), (kN|k1) [0228] the parameters Q, R, Q.sub.0.sup.1, R.sub.0.sup.1 (weighting matrices) and (regularization) [0229] 2. The calculation (k|k) is carried out, this corresponding to the estimation of the force applied by waves to the moving part, by performing the following steps: [0230] i. The following quadratic programming problem QP is resolved to obtain the solution *

[00042]$\{\begin{array}{c}\underset{}{\min }.Math.\left\{{}^T\left(P+{}_s^T.Math.R.Math..Math.{}_s\right).Math.-2\left[{\hat{w}}_0^T.Math..Math.{\hat{x}}_0^T.Math..Math.V^T\right]\left[\begin{array}{c}S\\ R.Math..Math.{}_s\end{array}\right].Math.\right\}\\ \mathrm{such}.Math..Math.\mathrm{as}.Math..Math.H_{}.Math.g_{}+H_d.Math.U_d\end{array}$ [0231] ii. the excitation force applied by waves incident upon the moving part w(k) is estimated by means of an equation with the form:

(k|k)=*(N).

[0232] 6) Controlling of the Wave-Energy System

[0233] During this step the wave-energy system is controlled so as to take into account the excitation force applied by the incident wave. This it is possible to drive the wave-energy system so as to optimize the recovered energy.

[0234] Controlling may consist of control of the moving part of the wave-energy system, for example by means of an electric, pneumatic or hydraulic conversion machine, called PTO (power take-off system). This PTO system influences the movement of the moving part and allows the mechanical energy to be transferred to the electrical, pneumatic or hydraulic network. Model predictive control (MPC) is an example of a method for controlling wave-energy systems.

Illustrative Example

[0235] In order to illustrate the advantages of the method according to the invention, the method according to the invention was tested in a tank. The wave-energy system corresponds to the wave-energy system shown in FIG. 1. The wave-energy system is equipped with a torque sensor for measuring the torque applied by the wave to the moving part. Beforehand, the torque applied by the wave is determined in the case where control is not active (namely the conversion machine does not generate any torque on the system) and the float is kept in an equilibrium position. These characteristics allows the excitation torque applied by the wave incident upon the moving part for this first test to be compared with other tests when the float is moving and the electric machine is activated.

[0236] For the present analysis a series of four types of different waves is considered. These four types of wave W1, W2, W3, W4 (from smallest to largest amplitude) are shown in FIG. 3 which illustrates the amplitude A as a function of the time T.

[0237] FIGS. 4a to 4c show a comparison of the excitation reference force (in this case the torque) applied by the wave incident upon the moving part, with the excitation force (in this case the torque) applied by the wave incident upon the moving part determined according to an embodiment of the invention. FIGS. 4a to 4c correspond to a wave of the type W1 shown in FIG. 3. The reference wave is indicated by REF and shown as a broken line. The forces of the incident wave determined by the method according to the invention are shown as a continuous line. FIG. 4a corresponds to the first embodiment of the invention MR1, FIG. 4b corresponds to the second embodiment of the invention MR2 and FIG. 4c corresponds to the third embodiment of the invention MR3. It can noted in these FIGS. 4a to 4c that, in all three embodiments, the graph representing the force determined by the method according to the invention is more or less superimposed on the graph representing the reference force. Thus, with the method according to the invention, it is possible to determine precisely the excitation force of the wave incident upon the moving part.

[0238] FIGS. 5a to 5c show a comparison of the excitation reference force (in this case the torque) applied by waves incident upon the moving part, with the excitation force (in this case the torque) applied by waves incident upon the moving part determined according to an embodiment of the invention. FIGS. 5a to 5c correspond to a wave of the type W3 shown in FIG. 3. The reference wave is indicated by REF and shown as a broken line. The forces of the wave determined by the method according to the invention are shown as a continuous line. FIG. 5a corresponds to the first embodiment of the invention MR1, FIG. 5b corresponds to the second embodiment of the invention MR2 and FIG. 5c corresponds to the third embodiment of the invention MR3. It can noted in these FIGS. 5a to 5c that, in all three embodiments, the graph representing the force determined by the method according to the invention is more or less superimposed on the graph representing the reference force. Thus, with the method according to the invention, it is possible to determine precisely the excitation force of the wave incident upon the moving part. It is also possible to note that the error between the reference force and the force which is determined is smaller for these graphs than for the graphs shown in FIGS. 4a to 4c.